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1.  Introduction 

Previously,  we  developed  a  method  to  allow  mammographic  differential  diagnosis  based 
upon  the  3-D  orientation  and  morphology  of  breast  calcifications  .  This  method  used  a 
limited- view,  binary  reconstruction  technique.  In  clinical  trials,  it  was  shown  to  be  of 
value  in  instances  where  calcifications  are  associated  with  a  mass.  In  such  cases,  we 
could  distinguish  between  preferentially  peripherally  distributed  calcifications  that  are 
predominantly  benign  and  homogeneously  distributed  calcifications  that  are  more  likely 
to  be  malignant.  We  have  also  been  able  to  elucidate  the  linear  distribution  of 
calcifications  contained  within  a  ductal  system.  Unfortunately,  this  reconstruction  method 
does  not  allow  one  to  image  non-calcified  tissues,  or  relate  calcifications  to  the 
surrounding  tissue.  Thus,  in  the  present  study  we  have  investigated  alternative  methods 
of  generating  3-D  images  of  the  breast,  namely  stereoscopy,  linear  tomosynthesis,  and 
micro-CT.  Such  methods  capable  of  producing  3-D  images  of  both  calcifications  and  the 
surrounding  breast  tissue. 

It  is  our  purpose  to  develop  a  method  that  will  be  clinically  viable  in  terms  of  dose,  image 
quality  and  equipment  cost.  We  believe  that  these  proposed  developments  will  further 
enhance  the  3-D  imaging  and  evaluation  of  breast  cancer  by  allowing  the  radiologist  to 
view  calcifications  in  relation  to  the  surrounding  tissue,  and  to  allow  3-D  imaging  of  non- 
calcified  breast  tissues  at  doses  which  are  clinically  acceptable.  Stereoscopy  has  the  value 
of  providing  depth  perception  of  tissues  with  little  additional  dose,  however,  often  the 
small  angle  separating  the  views  is  insufficient  to  completely  determine  the  causes  of 
superposition.  Tomosynthesis  requires  more  views  and  potentially  higher  dose,  but 
provides  better  separation  of  tissues.  Artifacts  from  the  reconstruction  algorithms  can 
blur  synthetic  tomograms.  CT,  while  providing  the  best  3-D  images  requires  doses  that 
are  not  clinically  acceptable. 

This  report  summarizes  the  work  on  this  grant. 
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2.  Body 

2.1 .  Summary  of  Work  Items 

It  is  useful  to  restate  the  work  items  listed  in  the  original  grant.  They  are  as  follows: 

Task  1 : Develop  3-D  imaging  techniques  (Months  1-24) 

Subtask  la:  Stereoscopy  (Months  1-4) 

Subtask  lb:  Linear  Tomosynthesis  (Months  5-8) 

Subtask  lc:  Limited  View  Reconstruction  (Months  9-20) 

Subtask  Id:  Computed  Tomography  (Months  18-24) 

Task  2:  Phantom  Development  (Months  1-24) 

Task  3 Evaluate  3-D  imaging  methods  with  phantoms  (Months  1-36) 

Task  4:  Acquire  image  datasets  (Months  13-36) 

Task  5:  Observer  study  (Months  25-36) 

2.2.  Development  and  Evaluation  of  3-D  Imaging  Techniques 

There  has  been  considerable  interest  in  imaging  the  breast  in  3-D.  This  has  included 
numerous  radiographic  methods  including  stereoscopy1’2,  tomosynthesis  ,  limited-view 
reconstruction  of  calcifications4  limited- view  tomography  reconstruction,  and 
computed  tomography7.  There  has  also  been  interest  in  3-D  ultrasound,  3-D  MRI,  and 
other  potential  methods  of  imaging  the  breast  in  3-D.  However,  to  date,  there  has  been  a 
sizeable  gap  between  the  proposed  techniques  and  clinical  feasibility.  To  achieve  clinical 
feasibility,  it  is  necessary  to  consider  these  techniques  on  the  basis  of  dose,  fundamental 
imaging  physics,  and  technology.  This  grant  was  designed  to  consider  4  such  methods 
and  compare  them,  both  on  the  basis  of  physical  performance  and  clinical  images. 

We  had  a  number  of  difficulties  completing  this  grant.  Significant  changes  in  clinical 
staff,  clinical  demands  on  the  time  of  the  participants  of  this  grant,  and  the  sale  of  a  key 
piece  of  equipment  by  Thomas  Jefferson  University  Hospital  all  served  to  retard  progress 
on  this  grant,  as  reported  previously.  In  spite  of  these  difficulties,  we  believe  that  we 
have  considerable  advanced  3-D  imaging  of  the  breast.  The  following  report  details  the 
results  of  this  grant. 


Image  Acquisition  Hardware 

The  grant  originally  proposed  to  use  two  pieces  of  hardware  for  the  completion  of  this 
grant.  The  first  was  a  prototype  device  built  in  Dr.  Maidment  s  laboratory.  It  consisted 
of  an  Eikonix  1412  linear  digital  camera,  coupled  to  an  x-ray  image  intensifier  (XRII)  via 
lenses.  These  were  mounted  on  an  optical  bench  with  a  Siemens  Bil50/30/50R  x-ray 
tube,  collimators,  and  a  specimen  holder  mounted  on  a  rotary  motion  stage.  A  Siemens 
Heliophos  5S  generator  supplied  the  x-ray  tube.  All  of  the  components  were  connected 
to  a  486-computer  running  Linux.  A  schematic  of  the  imaging  system  is  shown  in 
figure  1.  This  system  was  intended  for  preliminary  investigations.  It  had  known  image 
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Figure  1  A  schematic  of  the  imaging  system.  The  system  consists  of  an  x-ray  tube,  a  system  of 
translation  and  rotation  stages  for  positioning  the  specimen,  an  x-ray  image  intensifier  (XRII),  lenses 
and  a  CCD  camera.  Pre-  and  post-object  collimation  is  not  shown. 


quality  problems  that  would  prevent  it’s  use  through  the  entire  study.  Fortunately,  a 
second  system  existed  at  that  time  which  had  significantly  better  image  quality. 

The  second  system  consisted  of  a  Fischer  Mammotest  stereotactic  core  biopsy  system  and 
a  Fischer  MammoVision  digital  x-ray  detector  (Fischer  Imaging,  Denver,  CO).  We  have 
had  a  long-standing  relationship  with  the  Fischer  R&D  group,  and  Fischer  had  for  many 
years  provided  technical  support  for  this  work.  The  Fischer  system  was  supposed  to  be 
modified  to  include  a  computer-controlled  rotator  stage.  This  would  allow  us  to  mount  a 
specimen  holder  on  this  stage  and  use  the  Fischer  system  to  produce  tomographic  images. 
This  system  was  being  developed,  in  part,  because  of  its  physical  location  in  the  Breast 
Imaging  Center,  which  is  in  the  same  building  as  our  outpatient  breast  surgery  center. 

We  had  planned  to  interview  women  prior  to  needle  localization,  including  obtaining 
informed  consent.  Then,  following  their  surgery,  the  surgical  specimen  would  be  rapidly 
imaged  on  the  modified  Fischer  table.  Unfortunately,  the  hospital  sold  this  unit,  August 
10th,  1999,  just  8  weeks  after  we  hired  a  post-doc  for  this  grant. 

Since  that  time,  we  have  made  a  concerted  effort  to  allow  this  grant  to  proceed.  We  have 
done  this  in  two  ways.  First,  we  have  altered  the  original  imaging  system  by  replacing 
the  image  intensifier  and  Eikonix  camera  with  a  4”  XRII  (Toshiba),  and  an  SMD  (SMD, 
Fort  Collins,  CO)  1M30  CCD  camera.  The  acquisition  computer  and  electronics  were 
also  updated.  We  have  characterized  their  performance,  and  used  the  system  for 
computed  tomography.  Secondly,  we  obtained  a  DRC  (DRC  is  a  subsidiary  of  Hologic, 
Wilmington,  DE)  amorphous  selenium  active  matrix  detector.  This  system  was  used  for 
tomosynthesis  and  stereoscopy.  In  the  following  sections,  the  performance  of  the  various 
imaging  systems  is  discussed.  Then,  a  summary  of  the  research  in  the  four  imaging 
modalities  is  presented,  concentrating  primarily  on  stereoscopy. 
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Figure  2  Reconstructed  images  produced  with  the  Eikonix  detector.  The  left  image  is  a  wire  phantom 
designed  to  test  rotational  artifacts.  Note  that  the  wire  is  reconstructed  with  different  levels  of  success 
(they  should  look  like  dots,  not  crosses).  Also  note  the  circular  artifacts  due  to  the  instability  of  the  XRII. 
The  image  on  the  right  is  of  a  fresh  chicken  thighbone  (8.8  mm  x  7.4  mm),  clearly  showing  a  ring  of 
cortical  bone,  and  details  within  the  medullary  cavity.  This  image  still  demonstrates  subtle  circular 
artifacts  that  limit  the  detection  of  low  contrast  structures  in  the  bone. 


Detector  Characterization 

In  each  of  the  designs  discussed  below,  a  Siemens  Heliophos  5S  x-ray  generator  and  a 
Siemens  Bil50/30/50R  x-ray  tube,  a  Parker  computer  controlled  rotary  stage,  and  a 
custom  specimen/phantom  holder  are  used.  These  form  the  basis  of  the  imaging  system 
being  built.  They  are  all  mounted  on  a  optical  breadboard  that  is  equipped  with 
appropriate  rails  and  mounting  hardware.  A  rotate  only  geometry  is  used  to  acquire 
images.  The  setup  allows  acquisition  of  either  tomographic  or  stereoscopic  images, 
depending  upon  the  acquisition  protocol  used,  by  varying  the  angle  of  acquisition  and  the 
dose  per  image.  Each  detector  is  capable  of  acquiring  2-D  images,  and  hence  3-D 
volume  reconstruction  is  also  possible. 

Eikonix  Detector 

The  first  detector  built  consisted  of  an  Eikonix  1412  digital  linear  camera  and  a  Siemens 
9”/6”  x-ray  image  intensifier  (XRII).  Because  this  was  a  linear  camera,  true  CT  images 
were  acquired  differently  than  the  other  methods.  In  the  limited  view  methods  (including 
stereoscopy),  2-D  images  were  acquired  by  scanning  the  detector  at  selected  angles. 
However,  due  to  time  constraints,  CT  images  were  made  in  two  different  ways.  Angular 
sub-sampling  allowed  us  to  acquire  up  to  200  images  for  3-D  reconstruction. 
Alternatively,  1-D  images  could  be  acquired  at  a  greater  number  of  angles,  but  only  a 
single  2-D  slice  of  the  object  would  be  reconstructed. 
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This  camera  was  used  to  develop  the  image  reconstruction  algorithms  used  for  the  CT 
images  (discussed  in  detail  below).  In  the  original  grant  application,  the  CT  images  that 
were  shown,  were  of  very  high  contrast  objects,  and  were  acquired  as  1-D  samples  to 
produce  a  single  2-D  slice.  A  conventional  filtered  back-projection  reconstruction 
algorithm  was  used.  Since  that  time,  we  have  altered  our  image  acquisition  code  to  allow 
2-D  images  to  be  acquired,  and  now  use  a  modified  Feldkemp  algorithm  to  generate  a 
stack  of  2-D  slices,  which  are  rendered  as  3-D  volume  data  sets. 

Use  of  the  Eikonix  camera  was  discontinued  due  to  circular  and  cross  artifacts  illustrated 
in  the  reconstructed  images  shown  in  figure  2.  In  August,  2000,  we  decided  to  change 
detector  designs. 

PRC  Detector 

At  the  time  of  our  decision  to  abandon  the  Eikonix  camera,  we  considered  two  possible 
designs.  The  first  was  based  upon  a  DRC  active  matrix  detector  that  we  had  in  our 
laboratory.  This  system  has  very  high  modulation  transfer  function  (MTF)  (see  figure  3), 
and  excellent  noise  power  spectra  (NPS)  and  detective  quantum  efficiency  (DQE).  These 
data  were  measured  in  our  laboratory.  In  support  of  this  work,  we  did  extensive 
modeling  and  experimentation.  This  work  was  reported  in  Medical  Physics  in  October 
2000,  a  reprint  (Appendix  1)  of  which  is  included  with  this  report .  The  system  has  been 
used  to  acquire  tomographic  images.  Examples  are  shown  in  figures  8  and  9  during  the 
discussion  of  computed  tomography 

The  great  strength  of  this  detector  is  the  image  quality.  We  continue  to  use  it  for  at  least 
some  of  experiments  to  test  individual  3-D  acquisition  methods.  However,  this  detector 
is  very  slow,  producing  one  image  every  40  to  50  seconds.  Thus,  even  using  angular- 
undersampling,  tomographic  images  produced  from  200  individual  views  take  more  than 
2  hours  to  acquire.  Thus,  we  had  to  look  for  a  different  detector  to  perform  our  CT 
research. 
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Ip/mrm 


Figure  3  MTF  of  the  DRC  detector,  compared  with  the  ideal  MTF  of  a 
sine  function.  Note  that  a  small  amount  of  edge  enhancement  occurs  due 
to  the  use  of  an  amorphous  selenium  photoconductor 


SMD  Detector 


The  second  new  detector  design  was  based  upon  a  CCD  camera  that  was  coupled  with 
lenses  to  a  phosphor  screen,  or  to  a  newer  4”  Toshiba  XRII.  This  design  allows  us  to  use 
a  very  fast  camera.  The  camera  we  chose  is  capable  of  imaging  at  30  frames  a  second. 
This  allows  tomographic  images  to  be  acquired  in  as  little  as  7  seconds  for  200 
projections  to  34  seconds  for  1000  projections.  However,  now  we  faced  limitations  of 
the  x-ray  generator,  which  extend  that  time  to  a  few  minutes.  The  problem  with  this 
design  is  that  it  is  only  ever  going  to  be  useful  as  a  laboratory  system.  The  use  of  lenses 
to  couple  light  from  a  phosphor  screen  to  a  CCD  camera  is  quite  inefficient.  Our 


Figure  4  A  portion  of  a  resolution  phantom, 
showing  a  limiting  resolution  of  1 1 .3  lp/mm, 
acquired  with  the  SMD  detector,  a  Lanex 
phosphor  screen,  and  85mm:55mm  relay  lenses 
giving  an  overall  field  of  view  of  20  mm.  A 
subsection  of  the  image  is  enhanced  for 
printing. 
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calculations  have  shown  that  as  little  as  2  light  quanta  per  x-ray  interaction  may  be 
recorded.  We  have  shown  previously9'10  that  this  inefficiency  can  result  in  needing 
higher  doses  for  the  patients  (or  test  objects).  As  a  result,  we  could  only  image  tissue 
samples  and  inanimate  test  objects.  It  is  our  expectation  that  a  clinical  system  could  use 
one  of  the  newer  active  matrix  arrays  that  are  capable  of  10  to  30  frames  per  second,  and 
which  would  be  dose  efficient. 

We  obtained  a  demonstration  detector  from  the  manufacturer  during  the  period  of 
December  2000  to  February  2001,  and  finally  installed  the  permanent  system  in  May 
2001.  During  that  time,  we  wrote  the  necessary  software  to  control  the  camera  with  the 
existing  x-ray  system.  We  performed  preliminary  testing  of  the  camera,  and  were 
satisfied  with  the  performance.  An  example  of  a  resolution  pattern  is  shown  in  figure  . 

Image  Technique  Optimization 

Stereoscopy 

Stereoscopy  is  the  process  by  which  two  views  of  a  scene,  obtained  at  slightly  different 
angles  (such  as  that  due  to  the  displacement  of  our  two  eyes)  provides  the  viewer  with  the 
perception  of  depth  (i.e.,  the  ability  to  discern  that  one  object  is  behind  another).  This 
process  can  be  applied  to  radiographic  imaging  to  achieve  a  similar  perception  of  depth. 
Originally,  we  proposed  to  perform  a  Raleigh  discrimination  test  to  determine  optimal 
angular  separation  for  radiographic  imaging.  However,  work  by  others  in  the  field, 
including  Chan1  led  us  to  consider  other  pertinent  research  questions.  In  particular,  the 
issue  of  dose  was  particularly  important  and  previously  un-addressed.  We  have 
expended  considerable  time  on  this  question. 

In  the  case  of  a  quantum-noise  limited  detector,  signal  detection  theory  suggests  that 
stereoradiographic  images  can  be  acquired  with  one  half  of  the  per-image  dose  nee  e 
for  a  standard  radiographic  projection,  as  information  from  the  two  stereo  images  can  be 
combined.  Previously,  film-screen  stereoradiography  has  been  performed  using  the  same 
per-image  dose  as  in  projection  radiography,  i.e.,  doubling  the  total  dose.  In  the  resulting 
paper  (Appendix  5),  the  assumption  of  a  possible  decrease  in  dose  for  stereoradiography 
was  tested  by  a  series  of  contrast-detail  experiments,  using  phantom  images  acquired 
over  a  range  of  exposures.  The  number  of  visible  details,  the  effective  reduction  of  the 
dose  and  the  effective  decrease  in  the  threshold  SNR  were  determined  using  human 
observers  under  several  display  and  viewing  conditions.  These  results  were  averaged  over 
five  observers  and  compared  with  multiple  readings  by  a  single  observer  and  with  the 
results  of  an  additional  observer  with  limited  stereo  vision.  Experimental  results  show 
that  the  total  dose  needed  to  produce  a  stereoradiographic  image  pair  is  approximately  1.1 
times  the  dose  needed  for  a  single  projection  in  standard  radiography.  The  observed 
benefit  was  greatest  when  the  images  were  displayed  to  emphasize  the  quantum  noise, 
and  decreased  with  increasing  dose.  To  further  this  work  to  the  more  general  situation  o 
arbitrary  objects,  with  varying  degrees  of  depth  perception,  and  for  use  in  different 
imaging  systems,  it  was  necessary  to  develop  a  3-D  computer  breast  phantom.  This  is 
described  in  more  detail  below. 
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Linear  Tomosynthesis 


We  have  also  performed  preliminary  experiments  with  regard  to  linear  tomosynthesis, 
using  the  Fischer  Mammovision  stereotactic  breast  imaging  system,  prior  to  its  sale.  For 
each  tomosynthetic  data  set,  fifteen  images  were  acquired  as  the  x-ray  tube  was  moved 
through  a  30  degrees  arc.  The  advantage  of  synthetic  tomography  over  conventional 
tomography  is  that  the  set  of  15  images  can  be  used  to  reconstruct  an  arbitrary  number  of 
planes,  while  in  conventional  tomography  each  image  plane  requires  an  additional  x-ray 
exposure 

.The  principles  of  tomosynthesis  are  illustrated  in  figure  5.  As  the  x-ray  focus  is  moved, 
the  imaging  plane  P  is  held  fixed.  To  perform  a  reconstruction,  each  x-ray  image  is 
viewed  as  a  gray-scale  function  gi  defined  on  a  region  P f  of  the  imaging  plane  P ,  where  i 
enumerates  the  x-ray  exposure.  For  each  plane  Q  parallel  to  P  in  which  reconstruction 
will  be  performed,  each  projection  position  of  the  x-ray  focus  defines  a  one-to-one 
correspondence  of  points  in  the  plane  Q  with  the  points  in  the  plane  P.  Explicitly,  for  the 
z-th  position  of  the  x-ray  tube,  for  each  point  q  in  Q,  there  is  a  line  passing  though  the 
focus  and  q  which  meets  P  in  a  unique  point  p.  Thus  for  every  i  there  is  a  "projective" 
transformation  between  the  points  of  P  and  Q.  Further,  the  gray-scale  function  gi  on  Pr 
can  now  be  considered  as  a  function  g'i  defined  in  the  region  Qj  in  the  plane  Q,  and  the 
value  of  the  reconstructed  gray  scale  image  at  a  point  q  is  the  sum  of  all  the  g'i  which  are 
defined  at  that  point.  As  with  conventional  tomography,  for  objects  lying  within  the 
plane  Q,  the  functions  g'i  add  coherently  to  produce  a  focused  image,  while  for  objects 
outside  the  plane  are  blurred.  An  example  of  images  of  a  preliminary  tomosynthesis 
phantom  is  shown  in  figures  6  and  7.  The  phantom  consists  of  Lucite  spheres  and  cubes 
contained  in  a  water-filled  Lucite  box.  This  box  is  superimposed  upon  a  contrast-detail 
phantom.  In  figure  6,  three  projection  images  of  this  phantom  are  shown.  These  images 
demonstrate  the  overlaying  "clutter"  of  the  spheres  and  cubes  that  effectively  obscure 
most  of  the  features  of  the  contrast-detail  phantom.  Even  at  10  times  the  dose  this  image 
was  acquired  at,  we  could  not  visualize  more  than  5  objects  of  the  contrast-detail 
phantom.  In  figure  7,  we  show  three  reconstructions  of  the  phantom  at  different  heights 
above  the  detector.  Two  sections  contain  only  the  acrylic  spheres.  At  a  height  of  10 
mm,  however,  the  contrast-detail  phantom  can  be  seen.  In  this  image,  24  contrast-detail 
elements  are  visible.  The  increase  in  detection  is  due  to  the  reduction  in  the  overlaying 
clutter.  Thus,  reduction  of  this  structural  noise  results  in  an  increase  in  the  effective  SNR 
of  the  contrast-detail  elements,  without  requiring  an  increase  in  dose. 
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Figure  5  An  x-ray  tube  is  shown  at  three  locations  (1,2,3);  the  imaging  plane,  P,  is 
held  fixed.  An  image  of  the  letter  “A”,  located  in  a  plane  Q  is  projected  to  three 
different  locations  on  the  detector.  Note  that  the  individual  regions  Qi  overlap. 
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Figure  7  The  tomosynthetic  images  of  the  phantom  shown  in  figure  XX.  The  phantom  is  shown  at  three  different 
depths,  including  one  the  incorporates  the  contrast-detail  elements  (leftmost).  Numerous  elements  are  now  visib 


Computed  Tomography 

A  computed  microtomography  system  has  been  built  to  image  breast  specimens.  The 
system  has  evaluated.  As  discussed  above,  the  system  has  been  built  so  that  it  can 
accommodate  two  different  types  of  detectors,  a  DRC  active  matrix  detector,  and  one 
based  on  a  phosphor  screen  optically-coupled  to  a  CCD  camera.  The  DRC  detector  is 
complete,  while  the  CCD  based  detector  was  completed  in  the  summer  of  2002. 

We  evaluated  both  conventional  and  fully  3-D  CT  image  acquisition  methods. 
Experimentally,  we  have  acquired  images  using  between  200  and  1000  projections, 
rotating  over  an  angle  of  360°.  Given  the  acquisition  geometry  used,  there  is  a  minimum 
number  of  projections  that  are  required  to  avoid  undersampling  the  projection  space. 

This  number  depends  upon  many  factors  including  the  object  size  and  detector  pitch. 
However,  undersampling  generally  only  causes  artifacts  in  very  highly  attenuating 
objects  such  as  bone.  Thus,  given  the  low  attenuation  of  breast  tissues,  undersamphng  is 
likely  to  be  less  of  a  problem. 
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For  performing  simple  axial  reconstructions,  the  projection  data  are  reconstructed  using  a 
filtered  back-projection  algorithm  (RECLBL,  Donner  Laboratory).  A  Hanning  filter  is 
used  in  the  reconstruction.  Prior  to  reconstruction,  the  projection  data  are  corrected  for 
pixel-to-pixel  variations  in  the  detector  response,  fluctuations  in  x-ray  exposure  between 
projections,  and  error  in  the  center  of  rotation. 


Numerous  phantoms  have  been  constructed  and  imaged.  Two  are  shown  in  figure  8. 
Shown  are  a  tomographic  resolution  test  object  and  a  uniformity  test  object.  ^The 
resolution  test  object  has  5  rows  of  holes,  varying  in  size  from  1/16”  to  1/64”.  All  are 
clearly  visible.  The  uniformity  test  object  (right)  consists  of  a  Lucite  cylinder  1”  in 
diameter.  The  image  demonstrates  a  very  uniform  background  without  cupping  or  other 


Figure  8  Two  reconstructed  axial  images,  acquired  with  a  CT  scanner  based  on  the  DRC  detector  The 
image  on  the  left  shows  a  resolution  phantom,  with  elements  of  size  1/16”  to  1/64  all  visible.  On  the 
right  is  a  uniform  lucite  phantom,  clearly  showing  the  image  uniformity  and  low  noise  that  is  achievable. 


artifacts.  These  images  show  the  excellent  potential  of  using  the  DRC  detector  to 
produce  tomographic  images. 


Shown  in  figure  9  are  some  3-D  images  of  biological  tissue.  The  3-D  images  are 
rendered  in  two  different  ways,,  surface  rendering  and  volume  rendering.  Both  show 
exquisite  detail.  We  believe  that  these  are  the  first  CT  images  of periplaneta  Americana. 
We  chose  to  image  this  species  due  to  it  ready  availability  and  due  to  the  size  of  the 
anatomic  structures.  Simply  stated,  this  represents  a  hard  imaging  task.  Attention 
should  be  paid  to  the  quality  of  the  reconstruction  of  the  leg  muscles.  Note  that  the 
muscle  fibers  inside  the  insect’s  exoskeleton  are  clearly  visible  and  discemable. 
Similarly,  part  of  the  exoskeleton  is  seen.  The  muscles  are  less  than  1  mm  in  diameter 
and  less  than  200  microns  separate  them  at  their  closest  point. 
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Figure  9  Surface  rendered  (left),  and  volume  rendered  (right)  microtomography  images  of  species  periplaneta 
Americana  (common  American  cockroach),  acquired  with  the  DRC  detector.  The  insect  is  approximate  y 
mm  long.  The  muscle  fibers  and  the  exoskeleton  of  the  legs  are  clearly  visible. 


3-D  Breast  Model 

The  observer  experiments  for  stereomammography,  stimulated  considerable  interest  by  a 
number  of  groups,  and  motivated  us  to  consider  this  problem  more  carefully.  Much  o 
the  last  half  of  the  grant  was  spent  specifically  addressing  the  issues  of  the  dose 
requirements  for  3-D  imaging  of  the  breast.  Such  questions  must  be  addressed  either 
through  phantom  experiments  or  computer  simulations.  The  number  of  x-ray  images 
needed  to  perform  this  research  solely  using  phantom  images  was  prohibitive.  So,  we 
had  to  develop  a  suitable  computer  simulation  of  the  breast.  This  simulation  is  discussed 
in  detail  in  Appendices  2  and  3.  The  phantom  consists  of  a  3-dimensional  array  of 
geometric  objects  positioned  using  rules  developed  from  observing  MRI  and  subgross 
histology  specimens.  The  phantom  can  have  various  properties  associated  with  it,  such 
as  mechanical  properties,  x-ray  attenuation  coefficients,  MRI  or  ultrasound  properties. 

To  date  we  have  only  simulated  x-ray  imaging  with  breast  compression.  This  model  has 
generated  more  collaborative  research  interest  than  any  other  research  project  that  we 
have  ever  worked  on.  Currently  more  than  10  research  groups  are  using  this  model 
worldwide. 
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A  V 


Breast  Duct  Analysis 

Although  not  originally  part  of  the  statement  of  work,  there  has  been  an  interested  side- 
research  project,  which  developed  from  this  grant.  This  project,  based  upon  ramification 
matrices,  originated  when  we  found  it  necessary  to  simulate  the  breast  ductal  network  for 
the  breast  simulation  described  above.  We  evaluated  a  number  of  galactograms  to 
measure  the  branching  properties  of  the  ductal  tree.  These  measurements  were  then  used 
in  modeling  the  breast.  However,  it  was  noted  that  the  measured  ramification  matrix 
values  correlated  with  radiographic  appearance  of  the  galactograms.  We  have 
investigated  this  further,  and  have  shown  that  there  is  a  strong  correlation  between  ductal 
branching  and  disease  state.  This  work  has  been  published  (Appendix  4).  This  work  has 
let  to  more  collaborative  research,  were  we  intend  to  look  at  the  effect  of  breast  ducta 
network  structure  as  a  predictive  metric  of  breast  cancer  risk,  both  in  humans  and  munne 

models. 
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2.3.  Discussion  and  Summary  of  Scientific  Results 

This  grant  has  generated  significant  progress  toward  developing  an  experimental  and 
theoretical  (computational)  framework  for  performing  3-D  imaging  of  the  breast.  T  e 
experimental  framework  consists  of  an  x-ray  generator,  x-ray  tube,  collimators,  an 
rotary  stage  for  holding  specimens  or  phantoms.  Each  of  the  components  is  under 
computer  control.  The  framework  is  designed  to  use  standard  optical  mounting  hardware 
and  an  optical  breadboard  to  allow  testing  of  many  different  detectors  and  acquisition 

protocols. 

We  have  developed  and  characterized  two  detectors  that  fit  within  this  framework.  The 
first  detector  is  an  active  matrix  amorphous-selenium  device,  and  second  uses  a  phosphor 
screen  or  XRII  optically  coupled  to  a  CCD  camera.  The  first  device  has  excellent  dose 
efficiency  (DQE)  and  resolution,  but  is  too  slow  to  allow  easy  tomographic  imaging.  T  e 
latter  device  again  has  excellent  resolution,  but  lacks  the  dose  efficiency  of  the  flat  panel 
detector.  This  detector  does,  however,  have  the  temporal  response  needed  to  perform 

CT. 

We  have  performed  extensive  evaluations  of  the  reconstruction  methods.  In  particular, 
we  have  concentrated  on  stereomammography  and  addressed  the  issue  of  dose  m 
stereomammography.  We  have  shown  that  stereomammography  is  possible  at 
approximately  the  same  dose  as  conventional  projection  mammography  This  work 
considered  the  effect  of  x-ray  quantum  noise.  To  extend  this  work  to  other  3-D  imaging 
techniques,  and  to  consider  additional  effects,  such  as  visual  disparity  and  depth  ^ 
perception,  we  developed  an  extensive  3-D  computer  simulation  of  the  breast.  This  work 
in  turn,  has  led  to  efforts  in  characterizing  breast  tissue  based  upon  breast  ductal 

branching. 
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3.  Key  Research  Accomplishments 

The  following  is  a  list  of  key  research  accomplishments  resulting  from  this  work: 

•  Developed  a  generic  testbed  for  3-D  imaging  research,  consisting  of  x-ray  generator, 
rotary  stage,  specimen/phantom  holder,  and  detector  assembly. 

•  Developed  three  different  detectors 

o  the  original  was  based  on  an  Eikonix  linear  CCD  coupled  to  an  XRII 
o  the  second  was  based  on  a  DRC  active  matrix  detector 
o  the  third  was  based  on  an  SMD  30  fps  CCD  camera 

•  Developed  stereotactic,  tomosynthetic,  and  volume  rendering  display  software 

•  Developed  theoretical  framework  for  measuring  and  calculating  the  physical  performance 
of  aliasing  and  non-aliasing  digital  x-ray  detectors. 

•  Performed  dose  optimization  of  stereoscopy  based  on  experimental  and  theoretical 
techniques. 

•  Developed  a  3-D  breast  simulation  for  projection  and  tomographic  imaging  of  the  breast 
using  x-rays  (and  potentially  other  imaging  modalities). 

•  Developed  (as  a  result  of  the  model)  a  method  to  analysis  breast  ductal  network 
branching,  and  have  identified  a  relationship  to  radiographic  findings. 
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5.  Conclusions 

An  experimental  framework  for  performing  3-D  imaging  of  the  breast  has  been 
developed.  This  framework  consists  of  an  x-ray  generator,  x-ray  tube,  collimators,  and 
rotary  stage  for  holding  specimens  or  phantoms.  Each  of  the  components  is  under 
computer  control.  The  framework  is  designed  to  use  standard  optical  mounting  hardware 
and  an  optical  breadboard  to  allow  testing  of  many  different  detectors  and  acquisition 
protocols. 

We  have  developed  and  characterized  three  detectors.  The  first  is  now  obsolete.  The 
second  detector  is  an  active  matrix  amorphous-selenium  device,  and  third  uses  a 
phosphor  screen  or  an  XRII  optically  coupled  to  a  CCD  camera.  The  second  device  has 
excellent  dose  efficiency  (DQE)  and  resolution,  but  is  too  slow  to  allow  easy 
tomographic  imaging.  The  third  device  again  has  excellent  resolution,  but  lacks  the  dose 
efficiency  of  the  flat  panel  detector.  This  detector  does,  however,  have  the  temporal 
response  needed  to  perform  CT . 

We  have  performed  extensive  evaluations  of  the  reconstruction  methods.  In  particular, 
we  have  concentrated  on  stereomammography  and  addressed  the  issue  of  dose  in 
stereomammography.  We  have  shown  that  stereomammography  is  possible  at 
approximately  the  same  dose  as  conventional  projection  mammography.  This  work 
considered  the  effect  of  x-ray  quantum  noise.  To  extend  this  work  to  other  3-D  imaging 
techniques,  and  to  consider  additional  effects,  such  as  visual  disparity  and  depth 
perception,  we  developed  an  extensive  3-D  computer  simulation  of  the  breast.  This  work 
in  turn,  has  led  to  efforts  in  characterizing  breast  tissue  based  upon  breast  ductal 
branching. 
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The  optical  transfer  function  (OTF)  and  the  noise  power  or  Wiener  spectrum  are  defined  for 
detectors  consisting  of  a  lattice  of  discrete  elements  with  the  assumptions  of  linear  response, 
Gaussian  statistics,  and  stationarity  under  the  discrete  group  of  translations  which  leave  the  lattice 
fixed.  For  the  idealized  classification  task  of  determining  the  presence  or  absence  of  a  signal  under 
signal  known  exactly/background  known  exactly  (SKE/BKE)  conditions,  the  Wiener  spectrum,  the 
OTF,  along  with  an  analog  of  the  gray-scale  transfer  characteristic,  determine  the  signal-to-noise 
ratio  (SNR),  which  quantifies  the  ability  of  an  ideal  observer  to  perform  this  task.  While  this  result 
is  similar  to  the  established  result  for  continuous  detectors,  such  as  screen-film  systems,  the  theory 
of  discrete  lattices  of  detectors  must  take  into  account  the  fact  that  the  lattice  only  supports  a 
bounded  but  (in  the  limit  of  a  detector  of  arbitrarily  great  extent)  continuous  range  of  frequencies. 
Incident  signals  with  higher  spatial  frequencies  appear  in  the  data  at  lower  aliased  frequencies,  and 
there  are  pairs  of  signals  which  are  not  distinguishable  by  the  detector  (the  SNR  vanishes  for  the 
task  of  distinguishing  such  signals).  Further,  the  SNR  will  in  general  change  if  the  signal  is  spatially 
displaced  by  a  fraction  of  the  lattice  spacing,  although  this  change  will  be  small  for  objects  larger 
than  a  single  pixel.  Some  of  the  trade-offs  involved  in  detectors  of  this  sort,  particularly  in  dealing 
with  signal  frequencies  above  those  supported  by  the  lattice,  are  studied  in  a  simple  model. 
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I.  INTRODUCTION 

The  importance  of  signal  detection  theory  in  quantifying  the 
performance  of  medical  imaging  systems  (x-ray  screen-film 
imaging  being  perhaps  the  best  example)  gives  impetus  to 
applying  the  same  techniques  to  the  digital  radiographic  im¬ 
aging  systems  which  are  now  coming  into  clinical  use.  As 
applied  to  screen-film  systems,  signal  detection  theory  re¬ 
quires  three  assumptions  to  be  at  least  approximately  ful¬ 
filled:  that  the  detector  responds  linearly  to  the  incoming 
signal,  and  is  both  stationary  and  homogeneous  (i.e.,  both  the 
detector  response  and  the  additive  noise  are  translationally 
invariant).  One  can  then  summarize  the  response  of  the  sys¬ 
tem  in  terms  of  the  gray-scale  transfer  characteristic,  the  op¬ 
tical  transfer  function  (OTF),  and  the  noise  power  or  Wiener 
spectrum. 

The  digital  x-ray  imaging  systems  which  are  now  appear¬ 
ing  generally  behave  as  a  lattice  of  discrete  detector  ele¬ 
ments.  Although  digital,  these  detectors  are  generally  oper¬ 
ated  under  conditions  such  that  the  effects  of  quantization  are 
negligible.  When  compared  to  screen-film  systems,  these  de¬ 
tectors  tend  to  be  linear  over  a  wider  range  of  exposures. 
Like  screen-film,  for  low-contrast  signals  the  noise  is  ap¬ 
proximately  additive  and  Gaussian.  However,  as  the  size  of 
the  imaging  elements  is  now  comparable  to  the  size  of  some 
of  the  smaller  objects  which  are  of  clinical  interest  (around 
0.1  mm),  these  detectors  are  not  strictly  homogeneous  in  that 
translations  by  a  fraction  of  the  lattice  spacing  result  in  the 
signal  being  recorded  in  a  different  manner.  As  these  devices 
generally  consist  of  a  regular  lattice  of  sensitive  elements, 
they  still  possess  a  symmetry  with  respect  to  a  discrete  group 


of  translations.  This  symmetry  is  approximate  due  to  the  fi¬ 
nite  extent  of  physical  detectors.  However,  as  in  the  theory 
of  screen-film  systems,  corrections  for  the  limited  extent  of 
the  detector  are  negligible  for  many  practical  applications. 
Thus,  one  can  apply  Fourier  techniques  to  put  the  signal 
detection  theory  of  such  devices  in  a  form  which  is  both 
tractable  and  similar  to  the  theory  of  screen-film  systems. 
Instead  of  using  a  continuous  Fourier  transform,  one  uses  a 
discrete  space  Fourier  transform,  which  recodes  the  data  ac¬ 
quired  by  the  detector  at  a  discrete  lattice  of  positions  in 
terms  of  a  bounded  and  continuous  range  of  spatial  frequen¬ 
cies. 

For  screen-film  systems,  the  OTF  diagonalizes  the  linear 
operator  which  relates  the  input  signal  to  the  output.  As  de¬ 
tailed  below,  for  discrete- array  detectors  the  effects  of  alias¬ 
ing  introduce  a  null  space,  different  for  each  device,  which 
prevents  this  operator  from  being  diagonalized  using  a  basis 
common  to  all  devices,  but  the  OTF  represents  the  operator 
in  a  basis  in  which  it  is  sparse  in  the  sense  that  all  terms 
vanish  except  those  between  input  and  output  spatial  fre¬ 
quencies  which  are  equal  or  aliased.  The  Wiener  spectrum  is 
the  discrete  space  Fourier  transform  of  the  discrete  autoco¬ 
variance  function,  and  thus  is  also  defined  in  the  region  of 
frequency  space  which  the  lattice  supports.  As  in  the  case  of 
continuous  detectors,  for  low-contrast  objects  (so  that  re¬ 
sponses  are  approximately  linear),  these  quantities  determine 
the  signal-to-noise  ratio  (SNR)  which  is  an  appropriate 
figure-of-merit  for  the  classification  task  of  discriminating 
between  the  presence  or  absence  of  an  exactly  known  signal 
against  an  exactly  known  background  (SKE/BKE). 
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To  make  notation  definite,  it  is  necessary  to  review  the 
relevant  parts  of  signal  detection  theory  for  screen-film  im¬ 
aging  systems  (Sec.  II)  and  the  relevant  Fourier  techniques 
(Sec.  III).  Subsequently,  the  OTF  (Sec.  IV)  and  Wiener  spec¬ 
trum  (Sec.  V)  of  discrete  detectors  can  be  studied  leading  up 
to  the  calculation  of  SNR  for  the  case  of  SKE/BKE  (Sec. 
VI). 

As  discussed  below,  the  effects  of  aliasing  on  the  inter¬ 
pretation  of  the  OTF  have  been  noted  in  the  literature,  and 
the  definition  of  the  Wiener  spectrum  given  here  has  ap¬ 
peared  before.  However,  this  paper  presents  a  systematic 
theory  of  signal  detection  for  discrete-array  x-ray  detectors. 
The  approach  here  differs  from  the  more  general  theoretical 
approach  of  cross-talk  matrices'  in  that  the  work  reported 
here  uses  the  additional  assumptions  of  infinite  extent, 
Gaussian  statistics,  and  discrete  translational  symmetry. 
There  are  many  imaging  systems  for  which  these  assump¬ 
tions  are  not  appropriate  [e.g.,  three-dimensional  (3D)  to¬ 
mographic  reconstruction").  However,  the  work  presented 
here  has  certain  advantages  in  that  so  long  as  the  additional 
assumptions  are  approximately  true,  the  quantities  involved 
are  closely  related  to  those  used  with  screen-film  systems 
and  are  similar  or  identical  to  measured  quantities  relating  to 
digital  systems  which  have  appeared  in  the  literature. 

To  demonstrate  the  use  of  this  theoretical  structure  and  to 
investigate  some  of  the  trade-offs  inherent  in  discrete-array 
systems.  Sec.  VII  presents  a  simple  model  of  a  detector.  In 
particular,  this  gives  the  opportunity  to  investigate  how  the 
detector’s  response  to  high  spatial-frequency  signals  affects 
the  detection  of  objects  in  terms  of  SNR,  and  in  particular 
the  effect  which  is  sometimes  called  “noise  aliasing”  is  ad¬ 
dressed. 


II.  REVIEW 

The  theory  of  SKE/BKE  detection  and  classification  for 
screen-film-like  systems  has  been  expounded  in  detail,’ ~ 
and  is  reviewed  here  in  order  to  establish  notation  for  later 
comparison  with  discrete-array  detectors.  The  input  signal  is 
the  x-ray  intensity  per  unit  area  and  the  output  signal  is  the 
film  density,  both  as  a  function  of  position.  The  interesting 
and  tractable  case  is  the  search  for  low-contrast  variations 
/( r)  in  an  x-ray  beam  whose  baseline  intensity,  I„ ,  is  such 
that  the  changes  in  film  density  £>(r)  are  a  linear  function  of 
/( r).  The  additional  assumption  of  translational  invariance 
then  gives 

(D(r))=y(l0/gl—  f  fd2r'P( r-r')</(r')>,  (1) 

where  P(r-r'),  the  point  spread  function  (PSF),  is  the  den¬ 
sity  increase  of  the  film  at  position  r  due  to  a  given  x-ray 
intensity  at  r'.  The  brackets  (())  represent  ensemble  aver¬ 
ages,  i.e.,  averages  over  many  exposures.  Since  /(r)  and 
D(r)  are  here  defined  as  variations  relative  to  baseline  val¬ 
ues,  (/(r))  =  0  and  (D( r)}  =  0  in  the  absence  of  a  signal. 

Because  the  PSF  is  translationally  invariant,  the  convolu¬ 
tion  operator  is  diagonalized  in  frequency  space,  giving 
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(D(f))=——T——~  T(t)(I(f)),  (2) 

* O 

where  f  represents  a  two-dimensional  vector  in  frequency 
space  and  I ,  the  optical  transfer  function  (OTF),  is  the  Fou¬ 
rier  transform  of  the  PSF.  The  function  D( f)  is  the  Fourier 
transform  of  the  data  D(r),  following  the  convention  that  for 
any  suitable8  function  g, 

g( f)=  f  f  d2r  g(r)e~2viT  t, 


g(r)  =  J  J  d2f  g(f)e2mr  t,  (3) 

so  that  f  is  in  units  of  cycles  per  unit  length.  The  factor  of 
y(log|oe)//0  serves  to  convert  units  of  x-ray  beam  intensity 
into  units  of  change  of  film  density  in  the  region  of  linear 
response,  allowing  the  normalization  7(0)  =  1. 

Realizations  of  the  imaging  process  will  be  subject  to 
noise  which  can  be  characterized  by  the  autocovariance  func¬ 
tion 


C(r,  ,r2)  =  C(r2  ,r, )  =  (D(  r,  )D(r2)),  (4) 

that,  for  Gaussian  noise,  completely  determines  the  statis¬ 
tical  nature  of  the  noise  process.  For  stationary  processes,  the 
autocovariance  function  depends  only  upon  the  displace¬ 
ment, 

C(r1,r2)  =  C(r2-r1),  (5) 

so  that  the  autocovariance  C(r)  is  now  a  function  of  a  single 
vector,  the  displacement  r.  The  Wiener  spectrum  is  the  Fou¬ 
rier  transform  of  the  autocovariance  function,  and 


W({)=  lim  r~r 


J  J^d2rD(r)e 


- 2rrir •  f 


(6) 


shows  that  the  Wiener  spectrum  can  be  estimated  in  terms  of 
the  Fourier  components  of  signal-free  (“flat-field,,)  images 
over  regions  A  of  sufficiently  large  area  | A  | . 

For  Gaussian  noise  it  can  be  shown9  that  the  optimal  strat¬ 
egy  for  SKE/BKE  signal  detection  or  classification  consists 
of  choosing  a  mask  function  g(  r)  and  a  cutoff  value  for  the 
statistic 


6g  =  J  J  d2rg(r)D(r).  (7) 

The  efficacy  of  6g  for  discriminating  between  hypothesis  I 
(e.g.,  signal  absent)  and  hypothesis  II  (e.g.,  signal  present)  is 
measured  by  the  signal-to-noise  ratio 


SNR,; 


((<9g)i~(flg)n)2 
Var(  0g)  ’ 


(8) 


where  the  numerator  is  the  difference  between  the  expecta¬ 
tion  values  of  0g  under  the  two  hypotheses  and  the  denomi¬ 
nator  is  the  variance  in  the  statistic  6g ,  which  for  additive 
noise  is  independent  of  the  hypothesis. 

While  only  real- valued  functions  g(r)  are  needed  for 
calculating  decision  statistics,  it  will  be  useful  to  extend 
the  definition  of  6g  to  complex  valued  g(r),  in  which 
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case  the  variance  of  6g  is  the  sum  of  the  variances  of  the  real 
and  imaginary  parts.  The  variance  of  the  statistic  9g  is  given 
by 

Var(  0g)=((dg—(8g))($g—(6g))'*)  (9) 

=  J  f  d2rf  f  d2r'g(r)C(r,r')g(r')*  (10) 

=  J  f  d2fg(f)g*(f)W(f),  (11) 


as  can  be  seen  by  substitution  of  Eqs.  (7)  and  (4)  into  Eq.  (9) 
[indeed,  Eq.  (10)  can  be  taken  as  the  definition  of  the  auto¬ 
covariance  function].  It  can  be  shown,  for  example  using  the 
Schwarz  inequality,10'11  that  the  optimal  choice  of  mask 
function  is  given  by 


iz<*)= 


«/(f))n-</(f)>,mf) 

W(  f) 


(12) 


for  which  the  SNR  is  given  by 


y2(l°gio^)2  l(/(f))n-</(f)).l2lnf)|2 

I  a2  W(f) 


(13) 


where  the  subscript  J  indicates  that  this  represents  the  opti¬ 
mal  or  ideal7  observer  given  the  detector  and  task  at  hand. 

Returning  momentarily  to  the  task  of  estimating  the 
Wiener  spectrum  from  flat-field  images,  if  in  Eq.  (9)  one  sets 
g(r)  =  G(r)e2mt°'r>  where  G( r)  is  a  window  function  with 
normalization 


|  j  d2rG2(r)  =  l,  f  J  </2f|G(f)|2=  1,  (14) 

then  calculating  (Og0*)  by  Eq.  (11)  one  obtains 

f  f  </2f|G(f-f0)|2iy(f) 

=  ||  J  J  d2rG(r)e-2m'r°  rD(r)  j,  (15) 

which,  since  |G(f—  f0)\2  will  be  sharply  peaked  near  f0, 
shows  that  for  finite  length  data  sets  one  actually  estimates 
the  Wiener  spectrum  convolved  with  the  square  of  the  Fou¬ 
rier  transform  of  the  window  function.  In  particular,  for  a 
rect  window  so  that  G  is  chosen  to  vanish  outside  of  a  square 
region  of  area  |A|=L2  and  to  have  value  1/\J\A\  inside  that 
region, 

\rtt\P 

|c(f)|  -'H  ) 

= A  sinc2(L/v)sinc2(L/y),  (16) 

where  fx ,  fy  are  the  components  of  f.  For  large  areas,  this 
becomes  increasingly  like  a  delta-function,  giving  Eq.  (6). 


III.  MATHEMATICAL  PRELIMINARIES 

In  the  case  of  screen-film,  described  above,  we  studied  a 
mapping  of  functions  defined  for  all  spatial  positions  to  func- 
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tions  defined  for  all  frequencies,  namely  the  Fourier  trans¬ 
form  defined  by  Eq.  (3).  For  the  digital  x-ray  detectors  of 
interest  in  this  paper,  the  input  is  still  a  continuous  distribu¬ 
tion  of  x-rays,  but  the  output  signal  consists  of  data  with 
values  assigned  only  at  a  discrete  set  of  lattice  points.  The 
size  of  the  array  is  assumed  to  be  sufficiently  large  so  that 
boundary  effects  play  no  part,  and  thus  is  treated  as  being  of 
infinite  extent.  To  analyze  these  data  it  is  appropriate  to  use 
a  Fourier  technique  based  on  the  discrete  translational  sym¬ 
metry  of  the  detectors.  Mathematically,  this  transform  is 
similar  to  the  Fourier  series,  but  with  the  roles  of  the  space 
and  frequency  domains  reversed.12  Following  Giger,13  this  is 
called  the  “discrete”  space  Fourier  transform  (DFT).  The 
“finite”  Fourier  transform  (FFT)14  is  the  Fourier  technique 
applied  to  finite  sequences  of  data  points  that  is  customarily 
implemented  using  an  algorithm  known  as  the  “fast  Fourier 
transform.”  For  practical  purposes,  one  always  deals  with 
finite  data  sets,  and  the  discrete  space  Fourier  transform  is  a 
limiting  case  of  the  finite  Fourier  transform.  As  the  number 
of  equally  spaced  data  points  used  in  calculating  a  finite 
Fourier  transform  increases,  the  discrete  set  of  frequencies 
calculated  fill  more  and  more  densely  a  bounded  region  of 
frequency  space,  so  that  in  the  limit  one  obtains  a  function  of 
a  continuous  range  of  frequencies.  The  function  obtained  by 
this  limiting  process  is  the  discrete  space  Fourier  transform 
of  the  spatial  data. 

To  deal  with  the  two-dimensional  arrangements  of  sensi¬ 
tive  elements  which  are  of  interest,  it  will  be  convenient  to 
introduce  the  ideas  of  vectors  generating  a  lattice  and  of  dual 
basis  vectors.  A  two-dimensional  lattice  of  points  can  be 
specified  by  vectors  v{  and  v2  such  that  every  point  in  the 
lattice  can  be  represented  as 

rmi,m2=mly{  +  m2\2,  (17) 

where  m  x  and  m2  are  integers.  The  vectors  and  v2  are  said 
to  generate  the  lattice,  but  the  choice  of  vectors  for  a  given 
lattice  is  not  unique.  For  the  common  case  of  a  square  grid 
the  choice  of  Vj  as  lying  along  the  x-axis  and  v2  as  lying 
along  the  v-axis  is  natural.  The  plane  containing  the  lattice 
can  be  tiled  in  such  a  way  that  each  tile  contains  a  total  of 
one  lattice  point.  Each  tile  is  then  called  a  “unit  cell.”  For 
the  case  of  a  square  grid  of  detectors  with  spacing  a ,  the 
most  natural  choice  of  a  unit  cell  would  be  the  square  cen¬ 
tered  at  the  coordinate  origin  extending  to  ±a!2  along  both 
axes.  A  small  region  of  a  plane  containing  a  square  lattice 
and  one  of  the  unit  cells  is  drawn  in  Fig.  1(a).  The  reciprocal 
vectors  denoted  by  Wj  and  w2  are  defined  by  the  require¬ 
ments 


(18) 


and  serve  as  a  basis  for  frequency  space  and  as  generators  of 
the  reciprocal  lattice.  In  the  case  of  the  rectangular  grid  men¬ 
tioned  above,  each  w/  would  be  parallel  to  v;  and  scaled 
appropriately  as  illustrated  in  Fig.  1(b),  which  also  shows  a 
unit  cell  of  the  reciprocal  lattice.  To  help  clarify  these  ideas, 
Fig.  1(c)  shows  a  hexagonal  lattice  with  two  choices  of  the 
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Fig.  1.  (a)  a  rectangular  lattice,  (b)  The  reciprocal  lattice  of  (a),  (c)  A 
hexagonal  lattice,  (d)  the  reciprocal  lattice  of  (c).  Note  that  (a),  (b),  (c),  and 
(d)  represent  a  finite  region  of  a  lattice  which  covers  the  entire  plane,  (e)  A 
3X3  finite  rectangular  lattice,  (f)  The  circles  represent  the  frequencies  used 
in  the  finite  Fourier  transform  of  (e).  For  comparison,  a  unit  cell  of  the  full 
reciprocal  lattice  is  shown.  See  Sec.  Ill  for  details. 


unit  cell  (parallelogram  or  hexagon)  and  Fig.  1(d)  shows 
the  reciprocal  lattice  (note  that  v{  is  perpendicular  to  w2 
and  v2  is  perpendicular  to  Wi).  The  area  | A |  =  | V!  X  v2|  of  a 
unit  cell  is  independent  of  the  choice  of  unit  cell,  since 
it  is  fixed  by  the  average  density  of  lattice  points  over  large 
regions.  The  area  of  the  unit  cell  of  the  reciprocal  lattice, 
|AT|  =  |wj  X  w2|,  is  inversely  proportional  to  \A\,  as  can  be 
seen  by 


\A\\K\  = 


det 


det 


det 


(Vl).v 

(V2)x 

(v,M 

(V2)J 

Vl-W, 

vrw2\ 

v2*  Wj 

v2-w2/ 

,  o\ 
o  i ! 

=  1, 

(W|).v 

(Wl)y 


(W2).v\ 

(W  2)y) 


09) 


(20) 


making  use  of  the  fact  that  the  determinant  of  a  product 
of  matrices  is  equal  to  the  product  of  the  determinants  and 
Eq.  (18). 

For  any  function  g(nti,m->)  of  the  lattice,  the  discrete 
space  Fourier  transform  is  defined  *  as 


J(f>  =  |A |  2  g(mi,m2)e  2ir,r"'l.»'2'r  (21) 

W]  ,m2 

for  all  spatial  frequencies  f.  This  definition  is  equivalent  to 
evaluating  the  --transform  on  the  unit  circle  in  the  complex 
plane.12,15  It  is  also  equivalent  to  the  Fourier  transform  of 
the  function  obtained  from  the  data  set  by  interpolation  with 


sine  functions  (e.g.,  Ref.  16,  p.  230).  Direct  calculation  from 
Eq.  (21)  gives 

g(f)=g(f+m1W|+/M2W2),  (22) 

which  shows  that  g  is  periodic  in  Fourier  space  for  displace¬ 
ments  in  the  dual  lattice  and  one  need  only  consider  values 
of  |  on  one  unit  cell  of  this  lattice.  Any  frequency  f  outside 
of  this  unit  cell  is  an  alias  of  a  frequency  f'  inside  the  cell, 
with  f— f'  in  the  reciprocal  lattice.  Viewed  another  way,  the 
reciprocal  lattice  divides  points  in  the  frequency  plane  into 
equivalence  classes  of  points,  two  points  being  equivalent  if 
and  only  if  they  are  separated  by  a  vector  in  the  reciprocal 
lattice.  Any  unit  cell  will  contain  exactly  one  point  from  each 
equivalence  class  (except  for  boundaries),  and  knowledge  of 
g  on  the  unit  cell  determines  g  on  the  entire  plane.  Alterna¬ 
tively,  one  can  consider  g  as  being  defined  on  the  topological 
“quotient  space,”  a  torus,  just  as  one  can  consider  a  function 
on  the  real  line  with  period  2t t  as  defined  on  the  unit  circle 
(Ref.  17,  p.  155). 

The  exponential  functions  in  the  discrete  Fourier  transfor¬ 
mation  satisfy  a  simple  orthogonality  condition 

J J '  d2te~ (23) 

where  K  is  the  region  corresponding  to  the  unit  cell  of  the 
reciprocal  lattice  in  the  frequency  plane  and  |AT|  is  the  area  of 
this  region,  thus  giving 

g(nii  ,m2)  =  J  J  d2i  g(  f)e27rir'r,ni'm2  (24) 

as  the  inverse  transform.  The  complex  exponentials  form  a 
complete  set  of  orthogonal  functions,  so  that  any  appropriate 
periodic  function  of  frequency  f  can  be  represented  in  terms 
of  them.  The  completeness  can  also  be  expressed  in  terms  of 
a  comb  function  as 

2  2  <5(f-f'-f,,,,2),  (25) 

m  j  ’^2 

where  the  equality  is  interpreted  in  terms  of  distributions  and 
the  sum  on  the  right-hand  side  is  over  the  frequencies  in  the 
reciprocal  lattice.8,16,17 

For  actual  finite  data  sets,  one  applies  the  finite  Fourier 
transformation.  The  discrete  space  Fourier  transformation 
can  be  interpreted  as  a  limit  of  the  finite  Fourier  transforma¬ 
tion  as  the  number  of  equally  spaced  points  in  the  data  set  is 
increased.  Specifically,  consider  a  bounded  subset  of  the 
{r«, ,«,}  such  as  of  the  form 

•^={r„|.„2l V,=s«,<v;  ,N2^n2<N£,  (26) 

for  which  the  finite  Fourier  transform  and  its  inverse  are 
given  by 


gd\,h)=  2  .  g{nx,n2)e 


S(ni’"2)=AA^r( 


(27) 

(28) 


where  AV—V/-V,.  The  points  in  the  Fourier  space  are 
given  by 
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iW2 


X={th.h\L^ll<L\,L1*h<L'1},  (30) 

where  the  L’s  are  chosen  so  that  L/  —  Z^—A//  .  The  re¬ 

ciprocal  relationship  [Eq.  (18)]  between  the  basis  vectors 
{v,}  and  the  dual  basis  vectors  {wj  gives 


h  h  \ 

AA^Wl+AA^W2J'(niV,  +  "2V2) 


/  j/2  |  /2/22 

AA^+  AA^’ 


(31) 


which,  along  with  choosing  N\  =  —Ni=N0f2  and  L,=  —  L\ 
=N0/2  for  N0  even,  produces  a  more  conventional  represen¬ 
tation  of  the  finite  Fourier  transform. 

As  the  number  of  data  points  AA^AW?  increases,  the 
spacing  between  the  frequencies  j7  decreases,  so  that  in 
the  limit  the  data  points  on  the  lattice  extend  across  the  entire 
plane  and  the  frequency  values  fill  a  unit  cell  of  the  recipro¬ 
cal  lattice.  The  finite  sum  in  the  FFT  [Eq.  (27)]  approximates 
(with  a  factor  of  |A|)  the  infinite  sum  in  the  DFT  [Eq.  (21)], 
and  for  the  inverse  transform  the  sum  in  Eq.  (28)  (with  the 
introduction  of  a  factor  of  |A||Ar|  =  1)  becomes 


8  fft(wi>w2) 


y  W 

tlips^NAN2 


(|a|^ffx(/1  ,/2)) 


which  approximates  the  integral  used  in  the  inversion  of  the 
discrete  space  Fourier  transform,  Eq.  (24).  To  illustrate  this 
concept.  Fig.  1(e)  shows  a  small  rectangular  lattice  (corre¬ 
sponding  to  Ni=  —  1,  Nf;  —2).  The  circles  in  Fig.  1(f)  repre¬ 
sent  the  corresponding  frequency  vectors  for  use  with  the 
finite  Fourier  transform.  The  box  shows  the  region  which 
would  correspond  to  a  unit  cell  of  the  reciprocal  lattice  if  the 
lattice  in  Fig.  1(e)  were  extended  to  an  infinite  lattice.  If  the 
finite  lattice  shown  in  Fig.  1(e)  were  extended  (but  still  fi¬ 
nite),  the  corresponding  frequency  vectors  of  the  finite  Fou¬ 
rier  transform  would  fill  the  unit  cell  more  and  more  densely. 

It  should  be  noted  that  if  AA^  or  AV2  are  even,  then  some 
of  the  frequencies  at  which  the  finite  Fourier  transform  is 
defined  [shown  as  the  circles  in  Fig.  1(f)]  would  lie  on  the 
boundary  of  the  unit  cell,  and  such  frequencies  would  have 
aliases  which  also  lie  on  the  boundary.  For  example,  in  the 
square  lattice  considered  in  Fig.  1(f),  if  one  of  the  frequen¬ 
cies  at  which  the  finite  Fourier  transform  is  defined  fell  on 
the  edge  of  the  unit  cell,  an  alias  of  that  frequency  would  lie 
on  the  opposite  edge,  and  a  frequency  on  any  comer  would 
be  aliased  with  all  of  the  other  comers.  In  certain  sums  over 
frequency  components,  such  as  Eq.  (28),  it  is  useful  to  adopt 
the  convention  that  such  sums  include  exactly  one  represen¬ 
tative  from  each  class  of  aliased  frequencies,  so  that  frequen¬ 
cies  falling  on  the  boundary  of  the  unit  cell  are  not  counted 
multiple  times.  If  one  uses  the  “quotient  space”  point  of 


view,  this  follows  automatically  as  the  aliases  correspond  to 
a  single  point  in  the  quotient  space.  Alternatively,  one  might 
weigh  each  frequency  by  a  factor  (1/2  for  frequencies  lying 
on  edges  and  1/4  for  comers)  so  that  each  class  of  aliased 
frequencies  has  a  total  weight  of  1  (similar  to  counting  frac¬ 
tional  atoms  when  reckoning  the  number  of  atoms  in  a  unit 
cell  of  a  crystal). 

The  results  pertaining  to  Fourier  transformations  and  dual 
lattices  which  are  reviewed  in  this  section  have  direct  gener¬ 
alizations  to  any  number  of  dimensions,  but  as  the  statement 
of  the  results  for  arbitrary  finite  dimension  would  be  nota- 
tionally  cumbersome,  only  the  two-dimensional  results  have 
been  explicitly  stated.  For  notational  convenience,  let  m  rep¬ 
resent  the  ordered  pair  m1,w2,  so  that  g(m)  =  g(m1  ,m2) 
and  rm=r,fl|  „2,  and  similarly  for  k,  e.g.,  fk=ftl.*2. 

IV.  TRANSFER  FUNCTION 

The  analog  of  the  optical  transfer  function,  which  relates 
the  response  of  the  detector  to  the  input  signal  in  frequency 
space,  can  now  be  defined.  The  input  signal  (I)  is  a  continu¬ 
ous  function  of  the  plane.  As  (/)  is  defined  relative  to  the 
“flat-field,”  it  is  reasonable  to  assume  that  (I)  has  compact 
support,  or  at  least  vanishes  sufficiently  quickly  at  infinity  to 
leave  the  quantities  considered  here  well  defined.  Thus  the 
Fourier  transform  (/)  is  a  continuous  function  of  the  entire 
frequency  plane.  The  data  D(rm)  are  well-defined  only  at  the 
discrete  lattice  points  rm,  so  that  the  discrete  space  Fourier 
transform  (D( f))  is  determined  by  its  values  in  one  unit  cell 
of  the  reciprocal  lattice.  Values  of  ( D )  outside  of  the  first 
unit  cell  are  determined  by  the  periodicity  relative  to  the 
reciprocal  lattice  and  contain  no  new  information.  For  spatial 
frequencies  inside  the  first  unit  cell,  the  detector  responds  at 
the  same  frequency  as  the  input  signal.  For  frequencies  out¬ 
side  of  the  first  unit  cell,  the  detector  responds  at  an  aliased 
frequency,  so  it  is  impossible  to  uniquely  determine  the  input 
signal  without  additional  information,  although  it  will  be  ar¬ 
gued  in  later  sections  that  for  reasonable  tasks  this  is  not  a 
significant  problem. 

Each  point  on  the  detector  grid  is  assumed  to  respond 
linearly  to  the  incident  signal,  so  that  the  analog  of  Eq.  (1)  is 

<Z>(rm))  =  r J  f  d2r'P(rm, r')(/(r')>,  (33) 

where  P ,  the  analog  of  the  point  spread  function,  represents 
the  response  of  the  detector  at  rm  to  x-ray  light  incident  at  r\ 
and  T  is  a  constant  for  converting  x-ray  intensity  into  digital 
values,  generally  chosen  so  that  the  integral  of  P  with  respect 
to  r'  is  unity.  With  a  discrete  detector,  one  no  longer  has  full 
translational  invariance,  but  there  remains  an  invariance  un¬ 
der  translations  which  take  lattice  points  to  lattice  points, 
assuming  that  each  pixel  is  identical  except  for  position. 
Thus  we  can  write 

F(rm,r')  =  P(rm-r'),  (34) 

to  indicate  that  the  response  of  a  detector  element  to  an  input 
signal  depends  upon  the  displacement  of  the  detector  ele- 
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ment  from  the  region  to  which  the  signal  is  applied,  but  not 
upon  the  absolute  position  of  the  detector  element  or  the 
signal,  from  which  it  follows1819  that 

<D(rm))  =  rJ  I  d2r'P(rm—r')(I(r'))  (35) 

for  each  position  rm  of  the  sensitive  elements  on  the  lattice. 
While  the  data  D(rm)  are  only  available  at  the  lattice  points, 
the  convolution  can  be  calculated  at  any  point,  so  that 

^(r)  =  r|  |  dVP(r-r')</(r')>,  (36) 

^f)  =  rP(f)(/(f)>,  (37) 

serves  as  a  definition  of  3>{r)  for  any  position  r.  Although 
<%)(r)  is  equal  to  the  data  (D( rm)>  at  the  lattice  points  where 
r=  rm ,  at  other  points  @Kr)  is  an  interpolation  which  will  not 
in  general  represent  a  physical  quantity,  although  it  is  some¬ 
times  useful  to  think  of  &Hr)  as  the  response  of  a  virtual 
sensitive  element  added  to  the  detector  at  position  r  in  such 
a  manner  as  to  not  perturb  or  be  perturbed  by  the  other  ele¬ 
ments.  The  discrete  space  Fourier  transform  on  (Z)(rm))  can 
now  be  calculated  using  (D(  rm))  =  5"(rm)  for  g  ( rm )  in  Eq. 
(21),  giving 

<D(f)HA|2 

m 

=  |A|2  f  f  d2{’ &(f')e2TriTm'(-t' 
m  J  J  K 


=  S  ^(f+4)=rS  r(f+fk)</(f+fk)>. 

fk  fk 


which  follows  from  expressing  §0  in  terms  of  its  Fourier 
transform  and  using  the  completeness  relationship  expressed 
in  Eq.  (25). 

Comparison  of  Eq.  (38)  with  its  screen-film  analog, 
Eq.  (2),  helps  to  clarify  the  interpretation  of  the  OTF,  T(f). 
The  spacings  in  the  discrete  lattice  introduce  new  length 
scales  which  occur  explicitly  in  the  summation  over  aliases. 
In  the  limit  of  a  very  finely  grained  lattice,  so  that  | A  |  — >0, 
the  spacing  of  the  reciprocal  lattice  points  gets  larger,  until 
only  the  one  unaliased  term  contributes  significantly  to  Eq. 
(38),  and  the  screen-film  case  is  recovered. 

When  frequencies  higher  than  those  supported  by  the  lat¬ 
tice  are  present  in  the  signal,  the  summation  in  Eq.  (38) 
introduces  “aliasing,”  that  is,  there  exist  multiple  spatial 
input  frequencies  whose  output  is  at  the  same  frequency  and 
are  thus  not  distinguishable.  For  example,  considering  a  one¬ 
dimensional  lattice  with  pixel-pitch  of  1  cm,  oscillations  at  a 
rate  of  0.5  cycles  per  cm  can  not  be  distinguished  from  os¬ 
cillations  at  a  rate  of  1.5  cycles  per  cm.  From  Eq.  (38),  two 
components  of  the  input  signal  generate  the  same  component 
of  the  output  signal  if  and  only  if  their  spatial  frequencies 
differ  by  an  element  f*  <fc,  of  the  reciprocal  lattice. 

More  generally  for  any  lattice  there  are  frequencies  f0 
such  that  -{„  is  an  alias  of  fc  (for  example,  if  f  and  -f  are 
on  opposite  boundaries  of  the  first  unit  cell  in  the  reciprocal 


lattice).  For  such  a  frequency  f0  [noting  that  —  — 

and  (/(f»  =  (/*(-f)>  for  real-valued  P(r)  and  (/( r»]  it  is 
possible  to  choose  the  phase  of  (I ( f0 ) )  so  that 

r(f0)(i(fj)+T(-f0)(i(-f0))=o,  (39) 

showing  by  Eq.  (38)  that  a  sinusoidal  signal  concentrated 
at  frequency  f0  and  displaced  by  an  appropriate  offset  rela¬ 
tive  to  the  lattice  (as  determined  by  the  phase  of  (/(f0))) 
would  be  indistinguishable  from  the  flat-field  signal.  Return¬ 
ing  to  the  simple  one-dimensional  model  of  pixels  spaced  at 
1  cm,  this  result  means  that  for  some  displacement  relative 
to  the  lattice  the  input  of  a  sinusoidal  wave  of  frequency 
1  cycle/cm  would  give  vanishing  output.  If  the  detector  ele¬ 
ments  were  assumed  to  integrate  over  1  cm  intervals,  then 
the  output  vanishes  for  all  relative  phases  of  the  sinusoidal 
input  wave  and  the  lattice.  If,  alternatively,  the  detectors 
integrated  over  only  0.5  cm  regions  but  still  were  spaced  at 
1.0  cm  intervals,  then  the  sinusoidal  wave  would  have  van¬ 
ishing  output  only  when  the  nodes  of  the  sinusoid  fell  upon 
the  centers  of  the  0.5  cm  sensitive  regions  of  the  detectors 
and  would  otherwise  change  each  digital  value  by  a  phase- 
dependent  offset  from  the  flat-field  value. 

The  optical  transfer  function  has  been  written  in  terms  of 
a  Fourier  transform  using  complex  exponentials.  Since 
complex-valued  exponential  inputs  are  not  readily  available, 
it  is  necessary  to  ask  how  T  can  be  experimentally  measured. 
In  principle,  phantoms  machined  to  produce  sinusoidal  pat¬ 
terns  of  x-ray  intensity  could  be  used,  and  by  repeated  mea¬ 
surements  with  different  offsets  one  could  separate  the  posi¬ 
tive  and  negative  frequency  components.  A  more  practical 
method  is  the  well-known  slanted  edge  technique,20-21  in 
which  images  are  acquired  under  flat-field  conditions  except 
that  one  half-plane  of  the  detector  is  shielded  so  as  not  to 
receive  any  input  signal.  The  detector  response  D  as  a  func¬ 
tion  of  distance  from  the  edge  is  referred  to  as  the  edge- 
spread  function  ESF,  which  can  be  differentiated22  to  give 
the  line  spread  function,  LSF.  Alternatively,  by  providing  an 
appropriate  input  the  LSF  can  be  acquired  directly.2-  The 
LSF  represents  integrals  through  the  PSF  along  lines  parallel 
to  the  edge,  so  that  by  acquiring  data  with  the  edge  at  mul¬ 
tiple  angles  one  obtains  the  radon  transform  of  the  PSF.  One 
can  reconstruct  the  PSF,  but  it  is  more  common  to  stop  after 
computing  the  Fourier  transform  of  the  ESF,  which  gives 
values  of  the  OTF  for  spatial  frequencies  f  which  are  normal 
to  the  edge.  For  discrete-array  detectors  it  is  desirable  that 
the  slope  of  the  edge  is  not  commensurate  with  the  lattice 
spacing  (for  example,  on  a  square  grid,  if  the  edge  is  not 
parallel  to  one  of  the  axes  and  does  not  have  a  slope  which  is 
a  ratio  of  small  whole  numbers  like  1/2  or  2/3).  When  this 
condition  is  satisfied,  for  a  given  region  of  interest  the  dis¬ 
tances  z  of  the  lattice  points  rm  from  the  edge  will  be  distrib¬ 
uted  sufficiently  densely  and  evenly  so  that  the  ESF  is  said  to 
be  “super-sampled,”  i.e.,  sampled  at  a  rate  significantly 
higher  than  the  reciprocal  of  the  lattice  spacing,  so  that  it  is 
possible  to  measure  values  of  the  OTF  for  input  frequencies 
beyond  those  supported  by  the  lattice. 

For  discrete-array  detectors,  rotational  symmetry  will 
generally  be  only  approximately  valid  at  low  spatial  frequen- 
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ties,  so  it  is  desirable  to  make  measurements  at  multiple 
angles  relative  to  the  lattice.  When  the  ESF#  at  a  given  angle 
6  is  acquired,  it  is  often  the  case  that  the  precise  position  of 
the  edge  relative  to  the  lattice  is  not  known,  so  that  one 
actually  acquires  data  for  ESF^z  +  z#),  where  zq  represents 
the  lack  of  knowledge  of  the  exact  position  of  the  edge. 
Upon  taking  the  Fourier  transform  of  the  ESF,  this  intro¬ 
duces  a  phase  uncertainty  of  the  form  e2m^Zff  into  the  value 
of  r(f).  While  this  phase  uncertainty  also  occurs  in  measure¬ 
ments  of  screen-film  systems,  for  discrete-array  systems 
summations  over  aliased  frequencies  generally  are  not  pos¬ 
sible  given  uncertainties  in  the  relative  phases  of  values  of  T 
at  different  spatial  frequencies.  In  general  one  can  remove 
this  phase  uncertainty  by  redefining  the  lattice  positions  to 
correspond  to  the  “centers-of-mass”  of  the  response  func¬ 
tions  of  the  sensitive  elements.  More  specifically,  if 

J  J  d2rP(r)=J  dzESFff(z+zff)>0,  (40) 

then  it  is  possible  to  redefine  the  lattice  (by  a  shift)  so  that 
each  lattice  point  sits  at  the  center  of  mass  of  the  response 
function  of  the  associated  detector  element,  giving 

J  J  dxdy  xP(  r)  =  j  J  dxdy  yP{  r)  =  0.  (41) 

With  this  redefinition  of  the  lattice  position,  each  LSF  ac¬ 
quired  corresponds  to  a  radon  projection  of  the  PSF  (onto  a 
line  perpendicular  to  the  edge)  and  thus  the  center  of  mass  of 
each  LSF  should  be  at  the  origin.  This  corresponds  to  shift¬ 
ing  the  acquired  LSF  (adjusting  z$)  so  that 

J  dzESFe(z)z=0  (42) 

for  each  angle. 

As  a  practical  matter  this  results  in  an  increase  in  the 
amount  of  data  it  is  desirable  to  report  for  a  given  detector.  If 
one  can  assume  an  inversion  symmetry,  i.e.,  P(r)  =  P(  —  r), 
then  the  imaginary  part  of  the  transfer  function  will  vanish 
identically,  so  that  only  the  real  part  need  be  reported.  The 
absolute  value  of  the  OTF,  traditionally  called  the  modu¬ 
lation  transfer  function  (MTF),  gives  enough  information 
to  calculate  quantities  such  as  the  spatial  average  of  SNR2 
(Sec.  VI),  but  does  not  give  enough  information  to  explore 
other  aspects  of  the  detector,  such  as  the  spatial  variation  of 
SNR2  as  the  test  object  is  moved  relative  to  the  lattice.  Re¬ 
searchers  should  also  note  that  with  the  slanted  edge  tech¬ 
nique,  when  combining  raster  lines  to  plot  the  edge  spread 
function,  the  independent  variable  of  interest  is  the  distance 
from  the  edge,  which  for  square  lattices  differs  from  the 
distance  along  a  raster  line  by  a  factor  of  the  cosine  of  the 
angle  between  the  raster  line  and  the  normal  to  the  edge. 
This  factor  becomes  significant  when  trying  to  measure  the 
transfer  function  at  angles  away  from  the  detector  axes. 
Based  upon  the  experience  of  the  authors,  one  can  generally 
measure  values  of  the  OTF  at  frequencies  several  times  the 
highest  frequency  supported  by  the  lattice.  One  is,  of  course, 
measuring  the  response  of  the  detector  at  low  frequency 
aliases  to  higher  frequency  input  signals.  Whether  the  pres- 
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ence  of  these  aliased  signals  in  the  output  is  desirable  will 
depend  upon  the  task  at  hand.  For  example,  it  might  be  de¬ 
sirable  to  detect  a  high-frequency  signal  even  if  one  can’t 
distinguish  it  from  a  low-frequency  signal,  or  the  resulting 
ambiguity  might  be  unacceptable. 

The  question  of  what,  if  anything,  should  be  identified  as 
either  the  OTF  or  MTF  for  digital  systems  has  been  ad¬ 
dressed  in  several  ways  in  the  literature.  For  example, 
Dobbins24  discusses  the  “pre-sampled  OTF”  (OTFpre,  our 
T)  as  measured  via  the  LSF,23  but  then  emphasizes  the  fact 
that  the  response  to  an  input  signal  with  either  sinusoidal  or 
delta-function  spatial  variation  will  change  if  the  input  signal 
is  shifted  by  a  fraction  of  the  lattice  spacing.  This  depen¬ 
dence,  which  follows  from  Eq.  (38)  when  the  input  is  ex¬ 
panded  into  Fourier  components,  confounds  attempts  to  de¬ 
fine  the  MTF  either  in  terms  of  the  frequency  response  to  a 
single  delta  function  or  as  the  ratio  of  output-to-input  ampli¬ 
tude  for  a  sinusoid.  Dobbins  addresses  this  issue  by  defining 

OTF^(f)  =  2  OTFpre(f+  f|),  (43) 

f. 

and  defining  EMTF(f)  as  the  amplitude  of  the  detector  re¬ 
sponse  at  frequency  f  to  a  delta-function  input  averaged  over 
all  positions  of  the  delta  function.  Giger  and  Doi25  included 
such  a  summation  of  OTF  over  aliased  frequencies  in  their 
study  of  data  acquisition  and  display  for  digital  systems. 
Both  OTF^  and  EMTF  can  be  computed  in  terms  of  the  OTF, 
but  it  can  be  seen  that  neither  is  sufficient  for  calculating 
SNR.  Metz26  approaches  the  problem  in  essentially  the  same 
manner  as  discussed  in  this  paper,  and  indeed  Eqs.  (19)  and 
(30)  of  that  paper  essentially  give  our  Eq.  (38),  but  for  a 
slightly  more  specialized  case.  Metz  then  brings  up  the  point 
that  a  shift  by  a  fraction  of  the  lattice  spacing  in  the  input 
signal  does  not  result  in  a  simple  shift  in  the  output  data,  and 
concludes  that  “the  effect  is  accounted  for  mathematically, 
but  it  prevents  us  from  defining  a  unique  ‘transfer  function’ 
of  the  sampling  process. 

Experimentally,  Sones  and  Barnes27  recognized  the  desir¬ 
ability  of  measuring  the  transfer  function  above  the  maxi¬ 
mum  frequency  supported  by  the  sampling  lattice  in  their 
work  with  a  digital  radiography  unit.  This  measurement  was 
performed  using  a  novel  technique  based  upon  a  phantom 
consisting  of  periodically  arranged  wires,  the  distance  be¬ 
tween  the  wires  chosen  to  be  incommensurate  with  the  dis¬ 
tance  between  samples  acquired  by  the  detector.  Fujita,  Doi, 
and  Giger28  measured  the  “pre-sampling  analog  MTF” 
above  the  maximum  frequency  supported  by  their  sampling 
lattice  via  a  slanted  slit  technique  and  recognized  that 
“knowledge  of  the  pre-sampling  analog  MTF  ...  will  be  use¬ 
ful  in  the  determination  of  signal-to-noise  ratio  (SNR)  [and] 
the  evaluation  of  digital  systems,”  a  statement  with  which 
we  heartily  agree. 

Working  from  a  complementary  theoretical  perspective, 
Barrett  et alx  uses  the  “cross-talk”  matrix  to  address  the 
more  general  case  of  any  detector  whose  response  is  linear, 
then  proceeds  to  more  specialized  cases.  In  Barrett  et  al. ,  the 
input  to  the  system  is  defined  as  the  object  being  imaged 
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parameterized  in  terms  of  the  coefficients  of  its  three- 
dimensional  Fourier  series,  while  for  our  purposes  the  input 
is  the  x-ray  fluence  incident  on  the  detector.  For  projection 
radiography,  which  is  our  primary  interest,  the  incident  x-ray 
fluence  is  directly  related  to  the  integrated  attenuation  coef¬ 
ficient  of  the  object  along  rays  diverging  from  the  x-ray  fo¬ 
cus,  at  least  to  a  first  approximation.  As  our  goal  is  to  at¬ 
tempt  to  quantify  the  detector  response  independently  of 
other  technical  factors,  this  approximation  is  adequate.  Bar¬ 
rett  et  al.  is  concerned  with  detectors  which  may  have  rela¬ 
tively  few  sensitive  elements,  so  the  application  of  Fourier 
techniques  to  the  acquired  data  is  not  considered.  Barrett 
etal.  applies  the  cross-talk  matrix  to  the  case  of  a  one¬ 
dimensional  array  of  detector  elements  with  aperture  size 
equal  to  the  element  spacing,  and  finds  that  the  cross-talk 
between  components  of  the  input  at  separate  frequencies  de¬ 
creases  as  the  length  of  the  array  is  increased,  so  long  as  the 
frequencies  are  not  aliases  of  each  other.  Thus  in  the  limit  of 
a  homogeneous  detector  of  infinite  extent  one  recovers  the 
fact  that  the  transfer  function  behaves  as  a  sparse  matrix,  in 
which  all  terms  vanish  except  those  on  the  diagonal  or  relat¬ 
ing  aliased  frequencies. 

In  order  to  use  Eq.  (38)  to  calculate  the  response  of  the 
detector  to  a  given  input,  it  would  be  necessary  to  know  the 
position  of  the  object  being  imaged  with  a  precision  finer 
than  the  lattice  spacing.  Strictly  speaking,  to  calculate  the 
response  in  either  the  discrete  or  continuous  case  requires 
that  the  input  be  “perfectly  known.”  However,  in  the  case  of 
a  continuous  detector,  a  shift  in  position  of  the  input  will 
result  in  a  corresponding  shift  in  position  of  the  output,  while 
for  a  discrete  detector  the  “shape”  of  the  output  would 
change.  In  many  cases,  such  as  predicting  the  detectability  of 
randomly  placed  objects,  one  would  need  to  calculate  for  an 
ensemble  of  objects  displaced  with  random  phases  relative  to 
the  lattice,  as  will  be  illustrated  below  in  calculating  the  SNR 
of  small  objects. 


V.  NOISE 

Individual  realizations  of  an  imaging  process  have  an  ir¬ 
reducible  variability  which  sets  a  fundamental  limit  on  how 
effectively  the  detector  can  distinguish  between  various  in¬ 
puts.  For  discrete-array  systems,  as  for  screen-film  systems, 
the  noise  can  be  quantified  in  terms  of  the  autocovariance 
function.  If  the  noise  is  additive  and  Gaussian,  then  the  au¬ 
tocovariance  matrix  completely  summarizes  the  stochastic 
process  which  generates  the  noise.  If  the  system  is  also  sta¬ 
tionary,  then  Fourier  techniques  can  be  used  to  define  the 
Wiener  spectrum. 

The  discrete  autocovariance  function  is  given  by1  ' 

C(rm,r„)  =  <D(rm)D(rn)>,  (44) 

where  rm  and  rn  are  points  in  the  lattice  of  detectors,  the 
angled  brackets  represent  averaging  over  an  ensemble  of  flat- 
field  images,  and  as  discussed  above  (£>(rm))  =  0  in  the  ab¬ 
sence  of  a  signal.  Symmetry  under  interchange  of  positions 

C(rm,rn)  =  C(rn,rm)  (45) 
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is  an  immediate  result.  With  the  assumption  of  stationary, 
the  autocovariance  depends  only  upon  the  displacement 
rm — rn,  so  we  can  write 

C(rm,rn)  =  C(rm-r„)  (46) 

without  ambiguity.  Note  that  the  difference  between  two 
vectors  corresponding  to  lattice  points  is  again  a  vector  cor¬ 
responding  to  a  lattice  point,  so  C  on  the  right-hand  side  of 
Eq.  (46)  is  defined  at  precisely  the  lattice  points. 

The  Wiener  spectrum  W(f)  is  defined  as  the  discrete 
space  Fourier  transform  [Eq.  (21)]  of  the  autocovariance 
function  C(rm).  As  with  any  discrete  space  Fourier  trans¬ 
form,  the  Wiener  spectrum  is  periodic  in  frequency  space 
[Eq.  (22)]  so  that  one  need  only  consider  the  values  of  W(f) 
on  a  single  unit  cell  of  the  reciprocal  lattice.  It  is  noteworthy 
that  both  the  autocovariance  C  and  the  Wiener  spectrum  W 
are  real-valued  and  even.  As  with  screen-film  systems,  one 
considers  statistics  which  are  Unear  functions  of  the  data,  so 
if  g(ml,m2)  is  a  set  of  real  (or  complex)  numbers  defined  on 
the  lattice  points,  one  defines 

0g  =  lAl2  g(m)D(rJ.  (47) 

rm 

The  variance  of  0g  (for  g  complex  valued,  the  sum  of  the 
variances  of  the  real  and  complex  parts)  is  given  by 

Var(  dg)  =(0g0g)  (48> 

2  g(m)Z>(rm)j|2  g*(n)D(r „)j| 

(49) 

=  |A|22  2  g(m)C(rn-rm),g*(n)  (50) 

m  n 

in  terms  of  real  space.  Expressing  the  autocovariance  matrix 
as  the  inverse  discrete  space  Fourier  transform  [Eq.  (24)]  of 
the  Wiener  spectrum  one  obtains 

Var( 0g)=  |A|22  2  *(m>f  \/f 

m  n  J  J  K 

Xff(f)e2",(r,'r"y(n)  (51) 

=  f  jVf£(f)£*(f)W(f),  (52) 

where  the  second  step  follows  from  the  definition  of  the  dis¬ 
crete  space  Fourier  transform  [Eq.  (21)]. 

Thus  one  can  calculate  the  variance  of  a  statistic  0g , 
which  depends  in  a  linear  manner  upon  the  data,  using  either 
the  autocovariance  function  or  the  Wiener  spectrum.  Statis¬ 
tics  of  this  form,  for  g( m)  real-valued,  will  be  seen  to  cor¬ 
respond  to  decision  variables  of  ideal  observers  in  Sec.  VI. 
As  in  the  screen-film  case,  it  is  useful  to  consider  functions 
g( m)  corresponding  to  the  product  of  a  plane  wave  and  a 
windowing  function,  which  can  be  written  as 

gr  ( Ta)  —  G{m)e2m{°'Tm ,  (53) 

where  G(m)  is  a  real-valued  window  function  with  normal¬ 
ization 
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|A|2  G(m)G*(m)  =  l,  f  f  d2fG(f)G*(f)  =  1,  (54) 
m  J  J  K 

where  the  two  normalizations  are  equivalent  by  Parse val’s 
theorem.  Applying  Eq.  (49)  and  Eq.  (52), 


Var(0J  =  |A|2  2  G(m)e-2^o  ^D{rJ 


f  fKd2f\G(f-f0)j2W(f). 


For  suitable  windowing  functions  G,  |G(f— f0)|2  will  be 
strongly  peaked  near  f0  so  that  one  obtains  an  estimate  of  the 
Wiener  spectrum  at  the  specified  frequency,  W(f0).  In  par¬ 
ticular,  if  Grect(m)  is  chosen  as  \!{M  XM  2\A\)m  at  the  lattice 
points  me[0,...,Af 1]X[0,...,M2~  1],  then 

A/,- 1  M2- 1  j 

*»=.=  2  2  -f==f=fD( (57) 

mi=0  m2  =  0  \JM  2\ A  | 


|6(f-f*)l2= 


j£|  sin2(M,ir(f-f0)-v,) 

A/,M2  sin2(ir(f— fn)  •  v() 

sin2(M27r(f-f0)-v2) 

sin2(7r(f-f0)-v2) 


which  explicitly  shows  that  for  this  choice  of  G, 
\G(f-f0)\2  is  strongly  peaked  near  f0 .  For  a  square  lattice 
with  conventional  choice  of  basis  vectors,  (f —  f ,0)*v1  =  (/‘JC 
wher efx—(f0)x  *s  the  difference  in  the  x  com¬ 
ponents  of  the  frequencies  and  Ax  is  the  lattice  spacing 
in  the  x  direction,  and  similarly  for  the  y  axis.  In  general, 
if  a  separable  window  is  chosen,  so  that  G(m) 
=  GI(m1)G2(m2),  then  G(f)  =  G^f* v1)G2(f- v2),  so  that 
one  can  make  use  of  the  variety  of  one-dimensional  windows 
which  have  been  studied.29 

Returning  to  the  case  of  a  general  lattice,  Eq.  (58)  shows 
that  for  this  particular  choice  of  window,  as  is  typical,  the 
estimate  of  VP(f)  becomes  sharper  as  the  spatial  width  of  the 
window  increases,  so  that 


=  lim 

<Wtf,*2(  D), 

(59) 

M ,  M2- 

-vOO 

|A| 

m{m2 

A/,-1  M2~l  2 

2  2  D( rje-2’1'-*-  , 

m i=0  m2~0 

(60) 

where,  by  stationarity,  any  MxXM2  region  of  the  detector 
lattice  will  serve.  Specializing  to  the  zero  frequency  case, 
f=0,  one  gets 

W{  0)=  lim  (M{M2\A\) 

M  j  ,A/ 2~*co 

jj  {  A*,"  1  M2-\  \2\ 

j.D<r-,J  /•  <6,) 

which  is  the  discrete-array  version  of  Selwyn  granularity5 
(the  variance  in  the  average  digital  value  corresponds  to  the 


variance  in  the  spatially  averaged  optical  density  of  film). 
Comparing  Eq.  (61)  to  Eq.  (52),  one  can  interpret  Eq.  (61)  as 
the  statement  that  the  integrated  response  over  large  regions 
of  the  detector  depends  only  upon  the  low-frequency  com¬ 
ponents  of  the  Wiener  spectrum.  Viewed  spatially,  this  result 
means  that  the  digital  values  averaged  over  sufficiently  large 
disjoint  regions  are  approximately  independent,  so  that  the 
variance  of  the  average  over  N  large  subregions  scales  with 
l/N*l/MlM2. 

As  with  the  OTF,  the  results  of  the  screen-film  theory 
appear  as  a  limiting  case  for  sufficiently  fine  lattices.  Writing 
Eqs.  (59)  and  (60)  as 


W(  f)=  lim  •  -  |  j 

I  A/,-1  m2- i 

x(  2  2  \A\D{rm)e~2”itr« 

\  m |  =0  m2  —  0 


the  summations  become  approximations  of  the  integrals  in 
Eq.  (6). 

The  discrete  autocovariance  [Eq.  (44)],  the  definition  of 
the  Wiener  spectrum  as  the  discrete  space  Fourier  transform 
of  the  autocovariance,  and  the  use  of  Fourier  components  of 
flat-field  images  to  estimate  the  NPS  [Eq.  (59)]  have  oc¬ 
curred  in  several  places  in  the  medical  physics  literature,13  30 
but  historically  these  results  seem  to  have  been  considered 
less  than  satisfactory  from  a  theoretical  point  of  view.  For 
example,  Cunningham31  stated  that  while  “[i]t  is  tempting  to 
write  out  the  NPS  of  [the  sampled  digital  signal],  but  strictly 
speaking  this  violates  the  shift-invariance  assumption  since 
[the  data]  is  sampled  and  is  therefore  not  shift  invariant/' 
More  recently,  Cunningham,32  in  analyzing  the  concept  of 
NPS  in  terms  of  cyclostationary15,33  random  processes,  de¬ 
fines  WM^M^  [Eq.  (60)]  as  “a  working  definition  of  the  digi¬ 
tal  NPS."  As  detailed  in  Sec.  VI,  the  NPS,  as  defined  here,  is 
precisely  the  noise  which  sets  the  detection-theoretic  limits 
on  the  use  of  the  detector.  In  the  detection-theoretic  approach 
of  Barrett  et  a/.,1  the  Fisher  information  matrix  relates  the 
detector  noise  back  into  uncertainties  in  the  estimates  of  the 
Fourier  coefficients  of  the  object  being  imaged.  This  has  the 
advantage  that  it  removes  the  fundamentally  arbitrary  choice 
of  scale  in  using  digital  values,  but  if  aliased  frequencies 
become  important  the  Fisher  information  matrix  becomes 
singular  so  that  the  inversion  of  this  matrix  is  problematic. 

The  definition  of  NPS  given  here  is  intended  to  be  opera¬ 
tional  in  the  sense  that  it  is  defined  in  a  manner  which  can 
be  implemented  using  the  experimentally  available  digital 
values.  For  the  purposes  of  understanding  the  sources  of 
noise  in  detectors,  it  may  be  useful  to  consider  the  noise 
in  the  “presampled"  signal,  and  for  some  detectors  this  pre¬ 
sampled  signal  might  be  experimentally  accessible.  For  ex¬ 
ample,  in  a  detector  based  on  a  phosphor  screen  coupled  with 
a  lens  to  a  charge-coupled  device  (CCD)  camera,  one  could 
do  experiments  in  which  the  camera  is  replaced  by  a  photo¬ 
graphic  film.  For  some  devices,  such  as  TFT  arrays  using 
direct  conversion  mechanisms,  the  meaning  of  the  presa- 
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mpled  signal  is  less  clear  as  removal  or  refinement  of  the 
sampling  array  is  likely  to  change  the  electric  fields  respon¬ 
sible  for  charge  collection. 

As  reviewed  by  Wagner  and  Sandnk,30  the  calculation  of 
the  NPS  can  be  implemented  in  several  ways.  One  method  is 
to  estimate  the  autocovariance  function  [Eqs.  (44)  and  (46)] 
using  pairs  of  points  in  one  or  (preferably)  more  images,  and 
then  performing  the  Fourier  transform  to  give  the  NPS.  Al¬ 
ternatively,  the  variance  in  the  Fourier  components  is  used, 
as  in  Eq.  (55).  If  G  is  chosen  as  a  rectangular  window,  then 
Eq.  (55)  reduces  to  Eq.  (59),  so  that  <WM)A,2(f)>  [Eq.  (60)]  is 
used  as  an  estimate  of  W(f).  In  principle  the  frequency  f  is  a 
continuous  variable,  but  the  spread  of  |G(f)|“  limits  the  reso¬ 
lution  in  frequency  space  [by  Eq.  (56)]  and  this  spread  is 
inversely  proportional  to  the  size  of  the  spatial  region  and  on 
the  order  of  \K\!MxM2.  Given  this  resolution,  it  is  reason¬ 
able  to  calculate  the  NPS  at  M,M2  frequencies  spaced 
evenly  in  the  unit  cell  K  in  frequency  space.  Thus,  the  tech¬ 
niques  commonly  in  use  by  experimenters  give  precisely  the 
quantities  of  interest  from  our  current  theoretical  point  of 
view,  although  the  use  of  windows  other  than  the  rectangular 
window  might  be  of  interest  to  obtain  better  frequency  reso¬ 
lution. 

Generally,  frequency  resolution  is  not  a  limiting  factor  in 
estimating  the  Wiener  spectrum,  and  the  NPS  estimated  by 
( WM  w,(f)>  is  subjected  to  further  smoothing.  From  Eq.  (57) 
it  is  seen  that  is  the  variance  in  the  random 

variable  0rcct,  and  as  the  region  of  interest  used  in  the  cal¬ 
culation  is  made  larger,  the  variance  in  0lccl  tends  to  W(  f) 
which  will  be  nonzero  in  general.  Because  the  variance  of 
0rcc,  does  not  vanish,  neither  will  the  variance  in  |0recJ2,  so 
the  variance  in  VV'.w  f)  does  not  converge  to  zero  as 
Mx,M2— As  the  region  of  interest  is  made  larger,  one 
gains  in  spectral  resolution  but  not  precision,  and  this  repre¬ 
sents  an  unavoidable  trade-off.20'34  One  can  only  decrease 
the  uncertainty  in  the  estimates  of  the  Wiener  spectra  by 
averaging  estimates  of  W(  f)  from  several  different  regions  of 
interest.  Of  course,  for  the  purposes  of  analysis  one  could 
divide  a  large  region  into  several  smaller  regions,  and  the 
averaged  value  of  estimates  of  W(f)  would  then  have  less 
uncertainty,  but  the  spectral  blur  would  be  increased.  Since  it 
is  often  inconvenient  to  obtain  sufficiently  many  flat-field 
images  to  make  the  standard  error  in  the  estimates  of  VV(  f)  at 
individual  frequencies  small,  researchers  often  opt  for 
smoothing  the  experimental  spectrum. 

VI.  KNOWN  SIGNAL  DETECTION 

Having  addressed  the  issues  of  OTF  and  Wiener  spec¬ 
trum,  it  is  now  possible  to  use  the  signal-to-noise  ratio  (SNR) 
to  quantify  the  ability  of  the  detector  to  perform  SKE/BKE 
tasks.  First,  however,  it  is  useful  to  briefly  review  the  mean¬ 
ing  of  the  SNR  in  terms  of  an  ideal9'25  observer  working  with 
Gaussian  statistics.  The  ideal  observer  is  challenged  with  de¬ 
ciding  between  two  hypotheses  based  upon  a  given  set  of 
data.  In  the  current  context,  these  data  consist  of  the  digital 
values  obtained  from  the  detector,  and  for  the  moment 
we  will  restrict  the  observer  to  knowledge  of  only  a  finite 
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region  of  the  detector,  corresponding  to  indexes  meJ 
=  [MX,...,M ;-l]X[M2,...,M2-l].  This  observer  works 
under  the  assumption  that  given  hypothesis  Hx,  correspond¬ 
ing  to  an  expected  input  signal  (/(r))j  and  an  expected  data 
set  (£>(rm))i,  the  probability  density  function  describing  the 
expected  range  and  frequency  of  observed  data  sets  is  Gauss¬ 
ian.  This  Gaussian  distribution  in  (M[— Mi)X(M2-M2) 
=  AM,AM2  dimensions,  one  dimension  for  each  detector 
element  available  to  the  observer,  can  be  written  explicitly, 
but  to  make  the  formulas  somewhat  less  cumbersome  we 
use  the  following  notation:  Xm=/)(rm),  (■Xm)i=(£>(rm))|, 
(Xm)u=(D(rm))a,  and  {Xm}  =  {Z)(rm)|me..#}  is  a 
AM  |  A  Mi-dimensional  vector  in  the  space  of  all  possible 
data  values  for  the  detector  elements  in  region  yM.  The  prob¬ 
ability  distribution  which  governs  the  frequency  with  which 
particular  data  sets  will  be  obtained  under  hypothesis  Hx  is 
given  by 

P^S^X  ))=No£-1/22m.neVZ(Xm-(Xm)l)(C  ' *mn^n"(Xn)l)) 

ml'  (63) 

where  the  normalization  factor  is  given  by 


I  j  \(AM,4M2)/2  j 

N°  \2rrj  x/det  C  ’ 


(64) 


The  matrix  Cmn  is  the  autocovariance  function  C(rm,r„)  of 
Sec.  V  restricted  to  the  range  m,ne^.  The  fact  that  m  and 
n  are  double  indices,  e.g.,  m  stands  for  mx  ,m2,  is  not  a 
problem  from  the  theoretical  point  of  view,  and  in  principle 
for  a  numerical  calculation  one  could  simply  choose  a  con¬ 
venient  one-to-one  pairing  of  the  double  indices  mx,m2 
e./ffl  with  the  integers  1,...,AMiAM2  so  that  C  would  be 
indexed  in  a  more  customary  manner.  Under  hypothesis  H n, 
the  range  and  frequency  of  observed  data  sets  will  be  gov¬ 
erned  by  a  Gaussian  probability  density  P n ,  this  time  con¬ 
centrated  around  (X)a.  The  restricted  covariance  matrix,  C, 
occurring  in  both  cases,  will  be  the  same  under  the  assump¬ 
tion  that  the  noise  is  additive. 

Returning  to  the  question  of  how  to  decide  between  hy¬ 
pothesis  Hx  and  hypothesis  Hu ,  if  for  a  given  instance  of  the 
experiment  a  data  set  {Xm}={/9(rm)|me.  /X)  is  obtained 
such  that  Pu({Xm})  is  relatively  large  and  P[({Xm})  is  rela¬ 
tively  small,  it  would  generally  be  reasonable  to  favor  Hu. 
Thus  the  ideal  observer's  decision  rule  based  on  the  likeli¬ 
hood  ratio  Pn/Pi,  as  discussed  below,  is  intuitively  reason¬ 
able. 

The  ideal  observer  attempts  to  minimize  the  expected 
cost9,36  given  knowledge  of  the  cost  of  misclassification  un¬ 
der  either  hypothesis  and  the  a  priori  probabilities  associated 
with  each  hypothesis, 

(Cost)  =  P(Hi)  P(ChIl|l)  C^a+  P(Hn)  P(Chl|lI)Cn-i , 


where  in  the  first  term  P(  H j)  is  the  a  priori  probability  of  the 
state  corresponding  to  hypothesis  Hi  being  true,  P(ChIl|l)  is 
the  probability  of  mistakenly  choosing  hypothesis  Ha  when 
hypothesis  Hx  is  correct,  C^n  is  the  cost  associated  with  this 
error,  and  similarly  for  the  second  term.  Given  a  region  Rn 
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of  the  AM  i  AM 2  dimensional  data  space  and  the  decision 
rule  that,  if  the  observed  data  {D(rm)\me~/<i}  are  in  Ru 
then  the  observer  rules  in  favor  of  hypothesis  Hu  and  other¬ 
wise  in  favor  of  /7r,  then  the  probability  of  mistakenly  fa¬ 
voring  hypothesis  Hn  when  7/j  is  correct  is 

P(ChIl|l)  =  J  f  ...  d^^XjP^Xj),  (66) 

and,  as  under  either  hypothesis  the  total  probability  must  be 
unity, 

P(Chl|lI)  =  1  —  J  f  ...J*  dAM‘AMWPu({XJ)  (67) 

gives  the  probability  of  making  the  error  in  the  other  direc¬ 
tion.  Combining  Eqs.  (65)-(67), 

(Cost)  =  P(Hn)C  n_, 

+  f  f  fg  </“'“«DC({U),  (68) 


2  g„/Arm)  C(rm»rn)  =  ((77(ril))n— (£)(rn)),),  (74) 

meJ 

for  all  ney/^.  On  physical  grounds,  the  values  of  the  mask 
function  gJ6(  rm)  will  be  significant  only  in  the  region  near 
where  (7(r))n-(/(r))i  is  nonzero.  Further  from  this  region, 
the  values  of  gy^(r m)  will  tend  to  zero,  so  that  for  suffi¬ 
ciently  large  AMXAM2  the  ability  of  the  detector  to  dis¬ 
criminate  between  the  two  hypotheses  should  not  depend 
upon  the  exact  value  of  AM{  AM2  •  In  that  limit,  the  efficacy 
of  the  detector  for  the  SKE/BKE  task  should  be  set  by  the 
linear  statistic  6g  for  the  ideal  observer’s  mask  function  gx. 
This  mask  function  is  defined  implicitly  by 

\M  2  «2(rJC(rm,rn)  =  ((D(rn))n-(D(rn))I),  (75) 

me.// 

where  a  factor  of  |A|  is  introduced  to  simplify  the  form  of 
the  solution  which  in  the  Fourier  domain  is  given  by 


where 

DC({Xn})  =  P(Hl)Cl^Pl({Xj) 


<D(f)>n-<D(f)>i 

W(f) 


(76) 


-P(//„)CIMPn({Xm})  (69) 

is  the  differential  cost  which,  if  the  experiment  were  repeated 
sufficiently  often,  would  be  attributed  to  those  experiments 
which  gave  data  {D(rm)|me./^}.  Clearly  the  expected  cost 
given  by  Eq.  (68)  is  minimized  by  choosing  the  region  Rn  to 
be  precisely  the  region  where  the  differential  cost  DC  is 
negative,  so  that  the  ideal  observer’s  decision  rule  is  to 
choose  hypothesis  Hn  if  and  only  if  the  likelihood  ratio 

Pn({P(rJlnis^g}) 

P i({^(rm)|me-^})  '  ’ 

exceeds  the  threshold  value 


°  P(Hn)Cu. 


(71) 


By  adjusting  the  operating  point  A0  one  makes  trade-offs  in 
the  rates  of  the  two  possible  error  types,  as  can  be  shown 
graphically  in  terms  of  receiver  operator  curves  (ROC 
analysis).37  Equivalently  one  can  place  the  cutoff  on  log  A, 
and  from  Eq.  (63), 


log  A=  2  ((/>(rm))n-(7>(rm)),)(C-|)mn£>(rn) 

m,ne7/ 

+  const,  (72) 

where  the  constant  term  does  not  depend  upon  the  observed 
data.  Thus  an  ideal  observer,  viewing  a  finite  region  .  /6  of 
the  detector  array,  uses  a  linear  statistic  6//  defined  by 

Q.//-  2  g.//A*m)D(rm),  (73) 

m  g  „/6 

where  g  y  is  given  implicitly  by 


2f«/(f+  fk)>„- </(f+  fk)>,)  r(f+  fk) 

r — - - 

W(T) 


(77) 


The  statistic  0„  is  itself  a  Gaussian  variable  whose  variance 
can  be  computed  using  Eq.  (52),  so  that 


SNR|=r2 


f  ,  |2rk<A/(f+fk))r(f+fk)|2 

f/~(  W) 


(78) 


gives  the  SNR  corresponding  to  the  use  of  the  statistic,  as 
defined  in  Eq.  (8).  Thus  the  limiting  case  of  a  detector  array 
of  infinite  extent  is  well  defined,  for  pixels  “far  away’’  from 
the  region  of  interest  do  not  significantly  contribute  to  the 
decision.  Physically,  it  is  clear  that  the  “tails”  of  the  PSF 
and  autocovariance  functions  set  the  relevant  scale  by  which 
distance  from  the  edge  of  the  array  is  measured,  so  that  when 
the  projected  images  of  objects  appear  at  a  distance  from  the 
boundary  of  several  times  the  lengths  of  these  tails  the  de¬ 
tector  can  be  treated  as  essentially  infinite  and  Eq.  (78)  is 
valid. 

It  is  acknowledged  that  there  are  mathematical  subtleties 
related  to  a  truly  infinite  detector  which  are  not  addressed 
here.  For  example,38  the  data  set  for  such  a  detector  would 
represent  an  infinite  set  of  random  variables,  so  it  is  not 
possible  to  write  down  a  probability  density  distribution  like 
Eq.  (63)  in  the  infinite  case.  The  nature  of  the  physical  limit 
is  sufficiently  clear  that  a  study  of  these  mathematical  subtle¬ 
ties  could  not  change  the  results.  In  any  case,  the  fact  that  the 
linear  statistic  8g  with  g  =  gx  gives  the  optimal  SNR  of  any 
linear  statistic  can  be  proven  directly.  More  precisely,  if 
g(rm)  is  used  to  define  a  linear  statistic  6 g ,  then  letting 
AD  =  (D)U-(D)1, 
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J  J^d2fg(f)AD(f) 


J  JKd2f(g(f)sfim) 

J±m)2 

w W)l 


(79) 


where  the  second  step  is  an  application  of  the  Schwarz  in¬ 
equality.  Dividing  both  sides  of  Eq.  (79)  by  the  first  factor  on 
the  right,  one  obtains 


|//^2fg(f)AD(f)|2^  r  r  ,jAP(f)|2 
//^2f|g(f)|2lV(f)  JK  W(f) 


(80) 


where  the  quantity  on  the  left  is  the  SNR2  for  the  statistic  6g 
[Eqs.  (8)  and  (52)]  and  the  quantity  on  the  right,  proven  to  be 
larger,  is  the  SNR2  of  the  ideal  observer  as  given  by  Eq.  (78) 
[with  Eq.  (38)]. 

As  a  slightly  less  subtle  point,  the  construction  of  the 
ideal  observer  involves  dividing  by  W(f)>  which  is  problem¬ 
atic  if  W(f)  =  0  at  some  frequency.  For  physical  detectors, 
the  Wiener  spectrum  never  vanishes  as  there  is  always  some 
residual  noise.  Even  for  highly  idealized  detectors,  the 
Wiener  spectrum  must  reflect  the  noise  in  the  incident  x-ray 
fluence  so  that  it  can  only  disappear  at  frequencies  where  the 
OTF  vanishes,  and  at  these  frequencies  the  Wiener  spectrum 
will  vanish  no  faster  than  OTF2(f)  (discussed  in  more  detail 
in  the  next  section),  so  that  even  in  this  case  the  SNR  as 
given  by  Eq.  (78)  is  a  well-defined  limit. 

The  SNR  given  by  Eq.  (78)  corresponds  to  the  SKE/BKE 
decision  task  using  a  discrete-array  detector,  as  Eq.  (13) 
gives  the  SNR  for  the  SKE/BKE  decision  task  for  screen- 
film.  Strictly,  these  formulas  do  not  apply  to  the  task  of 
detection  when  the  observer  does  not  know  the  position  of 
the  object  being  imaged.  For  detecting  a  signal  of  unknown 
location,  one  can  calculate  the  ideal  observer's  SKE/BKE 
O^r)  for  each  possible  position  r  of  the  object.  A  common 
strategy  is  then  to  apply  a  threshold  to  Oj{r).  Under  the 
assumption  of  Gaussian  statistics  with  complete  knowledge 
except  for  position,  the  likelihood  ratio  computed  by  the 
ideal  observer  uses  0j(r)  in  a  nonlinear  manner35,29'40  that  is 
sensitive  to  peaks  in  0j{ r).  In  either  case,  the  values  of  SNR 
given  by  Eqs.  (13)  and  (78)  are  indicative  of  the  efficacy  of 
the  ideal  observer  in  the  more  general  case  of  the  position 
being  unknown. 

For  the  discrete-array  detector,  however,  the  value  of  the 
SNR  for  the  SKE/BKE  case  will  depend  upon  exactly  where 
the  object  is  relative  to  the  lattice.  While  this  variation  can  be 
significant  (for  example  detectors  could  have  interstitial 
spaces  where  objects  completely  disappear),  the  magnitude 
of  the  effect  decreases  for  objects  large  relative  to  the  lattice 
spacing.  Examples  of  this  for  several  simple  model  detectors 
will  be  given  in  the  next  section.  If  the  variation  in  SNR2 


with  position  is  not  too  great,  then  the  spatially  averaged 
value  of  SNR2  will  be  of  use.1  This  spatial  average  can  be 
computed  exactly  by  noting  that  if  an  object  is  shifted  by  a 
displacement  r,  the  Fourier  transform  is  multiplied  by  e2mtr 
so  that  in  Eq.  (78)  the  sum  over  elements  of  the  reciprocal 
lattice  becomes 

2 

2  (A/(f+f|c))r(f+f|c)e2mfk'r  ,  (81) 

fk 

where  a  common  factor  independent  of  k(|e27nf*r|  =  1)  has 
been  removed.  In  averaging  over  positions  r  in  Eq.  (78),  the 
denominator  of  the  integrand  does  not  depend  upon  r,  and 
the  numerator  is  the  square  of  the  magnitude  of  a  Fourier 
series  in  r,  so  that  in  integrating  over  r  to  obtain  the  average 
over  all  displacements  one  can  apply  Parseval’s  theorem  to 
obtain 


<SNR2)  =  r2J  j^d2 


2fjA/(f+fk)|2jr(f+fu)|2 

W) 


=r2/  /  <82) 

where  the  second  step  follows  from  noting  that  the  sum  of 
the  integrals  over  each  unit  cell  is  equivalent  to  the  integral 
over  the  entire  plane. 

As  for  the  OTF  and  NPS,  the  film-screen  result,  Eq.  (13), 
can  be  recovered  from  the  discrete-array  result  [Eq.  (78)]  by 
going  to  the  limit  of  a  sufficiently  fine  lattice,  in  which  case 
the  distance  to  the  first  aliased  frequency  is  so  large  that  only 
the  unaliased  term  contributes  to  Eq.  (78).  Similarly,  for  a 
sufficiently  fine  lattice  all  objects  are  large  relative  to  the 
lattice  spacing,  so  that  SNR2  does  not  vary  appreciably  as  the 
object  is  moved  relative  to  the  lattice  spacing.  These  facts 
prompt  the  identification1  of 

GNEQ(f)  =  T2|  r(f)  |2<f>2/  W(f) ,  (83) 

as  a  generalization  of  the  concept  of  noise  equivalent  quan¬ 
tum  flux  (NEQ),  where  <I>  is  the  incident  x-ray  flux,  and 

GDQE(  f )  =  T 2 1  T{  f)  | : 20>/  W(  f) ,  (84) 

as  a  generalized  detective  quantum  efficiency  (DQE).  These 
results  parallel  the  screen-film  theory,  except  that  factors  of 
fluence  appear  in  the  numerator  as  the  response  of  digital 
detectors  is  linear  with  fluence  [Eq.  (33)]  while  film  density 
is  linear  with  respect  to  the  log  of  fluence  [Eq.  (1)].  While 
Eq.  (82)  is  exact  in  the  context  of  the  assumptions  we  have 
made  about  the  detector,  SNR2  enters  nonlinearly  into  other 
quantities  such  as  the  various  probabilities  of  misclassifica- 
tion  for  a  given  operating  point  (sometimes  called  the  false 
positive  fraction  and  the  false  negative  fraction  in  ROC 
methodology).  However,  when  the  variation  in  SNR2  is  not 
too  large,  perhaps  as  measured  by  the  rms  (root-mean-square 
variation  in  SNR2),  then  the  GNEQ  and  spatially  averaged 
SNR2  can  be  considered  a  useful  summary  of  the  efficacy  of 
the  detector. 

In  this  paper  we  have  applied  the  concept  of  an  ideal 
observer  directly  to  the  digital  data.  The  results  obtained  are 
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implicit  in  the  work  of  Giger  et  a/.,13,25’41,42  but  Giger  et  al 
concentrates  on  issues  of  display  and  models  of  human  vi¬ 
sual  response  to  the  displayed  data.  As  these  tasks  are  decou¬ 
pled  from  image  acquisition  for  digital  systems,  it  is  worth 
considering  figures  of  merit  for  the  data  acquisition  system 
independent  of  the  display,  as  done  here.  The  results  of  this 
section  also  follow  as  limiting  cases  of  the  work  of  Barrett 
et  al.1  Of  particular  note.  Sec.  V  A1  discusses  a  simple  bin¬ 
ning  detector  and  obtains 


M 


SNR2 


A  fm 


m  =  1  CT 


(85) 


where  crm  is  the  uncorrelated  noise  in  the  mth  detector  and 
A gm  is  the  expected  change  in  the  data  value  at  the  mth 
detector  which  would  be  caused  by  the  signal.  This  particular 
result  can  be  obtained  directly  from  first  principles  based  on 
counting  statistics  in  each  detector  element.  In  the  stationary 
case,  (Tm  =  ar  is  a  constant,  so  in  Eq.  (78)  W(f)  =  a2\A\  and 
the  numerator  [using  Eq.  (38)  and  Parseval’s  identity]  be¬ 
comes  |A|2m|A(Z)(rm))|2,  again  recovering  the  result  [Eq. 
(85)]  based  on  counting  statistics  for  uncorrelated  noise.  It  is 
worth  noting  that  if  one  does  not  include  the  aliased  terms  in 
the  numerator  of  Eq.  (78)  (perhaps  on  the  grounds  that 
aliased  signals  are  not  useful),  the  value  of  SNR2  will  be 
underestimated.  The  aliased  response  is  part  of  the  physical 
response  of  the  detector,  and  in  this  case  the  aliased  terms 
will  add  coherently  in  such  a  manner  as  to  bring  the  calcu¬ 
lated  value  of  the  SNR2  up  to  the  value  in  Eq.  (85)  obtained 
from  counting  statistics. 


VII.  MODEL  DETECTORS 

To  give  a  feel  for  the  implications  of  the  above  theory,  the 
capabilities  of  detectors  with  reasonably  realistic  parameters 
will  now  be  investigated.  The  modeling  is  somewhat  simplis¬ 
tic,  but  sufficient  to  demonstrate  several  interesting  proper¬ 
ties,  such  as  the  dependence  of  SNR  on  the  position  of  the 
object  being  imaged,  and  certain  trade-offs  inherent  in  such 
detectors,  particularly  those  trade-offs  related  to  the  possible 
suppression  of  input  spatial  frequencies  above  the  frequen¬ 
cies  supported  by  the  lattice.  The  incident  x-ray  fluence  <5 
has  a  white  Wiener  spectrum,  W,(f)  =  <I>.  Among  other  sim¬ 
plifications,  which  will  be  discussed  in  more  detail  at  the  end 
of  the  section,  we  assume  100%  of  the  x-rays  interact.  Each 
x-ray  undergoes  a  stochastic  amplification,  characterized  by 
an  average  of  m  secondary  quanta  per  x-ray  with  am  =  yjm 
for  a  Poisson  process,  and  the  secondary  quanta  undergo  a 
stochastic  scattering  process,  with  a  spread  function  Ps  and 
transfer  function  Ts ,  before  being  “binned”  by  the  detector 
elements.  The  result  is  an  average  of  md>  secondary  quanta 
per  unit  area  on  the  detector  with  a  pre-sampled  Wiener 
spectrum  given  by43 

WJ(f)=[m2iy,(f)  +  ‘l>^-m<D]|rj(f)|2+/n<D.  (86) 

For  a  square  lattice  with  spacing  L,  binning  can  be  consid¬ 
ered  as  a  deterministic  convolution  with  rect  functions  rep¬ 
resenting  the  detector  regions,  so  that  with 


Tb(  f)  = 


sin(  ttL/x)\ 

-*Lfx  I 


sin(  TrLfy)\ 

*Lfx  r 


the  digital  noise  power  spectrum  can  be  written 


(87) 


W(t)=  As  |A|2a>(m2|7'J(f+fk)|2  +  m)|rfr(f+fk)|2 

m  fk 


+  ^£,  (88) 

where  the  factor  of  1/m2  is  introduced  so  that  digital  values 
will  correspond  to  x-ray  count  and  WE  is  the  electronic 
noise.  With  the  present  conventions  the  gray-scale  character¬ 
istic  is  set  to  T  =  |A|.  A  simplification  can  be  achieved44 
using 


sin(7r(x  +  «)) 
7 T(x+n) 


(89) 


for  any  x,  which  can  be  proven  by  applying  Parseval’s  theo¬ 
rem  to  the  Fourier  series  for  e2mxy  for  y  e  [  —  0.5, 0.5].  The 
experimentally  observable  transfer  function  (as  obtained,  for 
example,  by  the  slanted-edge  technique,  cf.  Sec.  IV)  contains 
the  effects  of  stochastic  scatter  and  binning,  thus  T{  f) 
=  Ts(t)Tb(f ),  so  that 

W(f)/|A|  =  <D|A|2  |r(f+fk)|2+^+HV|A|  (90) 

fk  m 

is  the  Wiener  spectrum  of  the  model  detector,  with  the  aver¬ 
age  number  of  x-rays  per  pixel  being  <3>|A|.  The  summation 
over  aliases  in  Eq.  (88)  is  often  referred  to  as  “noise  alias¬ 
ing.”  The  division  into  aliased  and  unaliased  components  is 
useful  for  modeling  a  variety  of  detectors,  but  it  should  be 
noted  that  this  division  is  generally  not  directly  experimen¬ 
tally  accessible,  at  least  not  without  modifying  the  detectors, 
and  that  in  principle  there  could  be  devices  which  are  sta¬ 
tionary,  and  therefore  have  Wiener  spectra,  but  for  which  the 
division  of  the  NPS  into  aliased  and  unaliased  components  is 
not  useful. 

It  is  useful  to  choose  values  of  the  parameters  in  the 
model  which  are  representative  of  detectors  of  current  clini¬ 
cal  interest,  as  this  can  help  in  the  understanding  of  the  phys¬ 
ics  which  determines  the  performance  of  these  devices,  but 
detailed  modeling  for  quantitative  comparison  to  actual  de¬ 
vices  is  beyond  the  scope  of  this  article.  We  assume  a  square 
lattice  with  spacing  of  L  =  0.143  mm,  operation  at  an  expo¬ 
sure  corresponding  to  |A  |<I>  =  1400  x-rays  per  pixel,  and  an 
amplification  factor  of  m  =  1000.  For  the  stochastic  transfer 
function  Ts  we  consider  three  possibilities:  a  “blur-free” 
detector  for  which  T5(f)  =  l,  typical  of  photoconductive 
arrays,45  and  two  “alias-free”  detectors  whose  stochastic 
transfer  functions  are  of  the  form  r?(f)  =  £~x^  with 
X  =0.463  and  \=0.34  mm,  which  approximates  the  transfer 
function  for  evaporated  Csl.46-48  Typically  electronic  noise 
}(We)/\A\  is  on  the  order  of  3-5  x-rays,  so  values  of  0,  42, 
and  8 2  cover  the  range  of  values  for  <F|A|/m  +  WEf\A\. 

The  transfer  functions  for  these  models  are  shown  in 
Fig.  2.  As  the  pixels  are  symmetric  with  respect  to  inversion 
through  their  centers  [i.e.,  for  the  PSF,  P(r)  =  P(-r),  and 
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Transfer  Functions 


GDQE 


Fig.  2.  The  optical  transfer  functions  of  three  model  detectors.  The  “blur- 
free  ’  '  detector  bins  the  secondary  quanta  without  smoothing,  while  for  the 
“alias-free"  detectors  the  distribution  of  secondary  quanta  is  smoothed  by 
an  exponential  MTF  (<TX||1)  before  binning.  Data  are  shown  as  a  function  of 
the  magnitude  of  the  spatial  frequency  for  several  angles. 

P(r)  is  a  real  number],  the  imaginary  part  of  the  transfer 
function  is  identically  zero,  so  only  the  real  part  need  be 
graphed.  The  OTF  is,  of  course,  a  function  of  two  variables, 
f  x  an dfy .  To  show  this,  we  plot  the  OTF  as  a  function  of  the 
magnitude  of  the  frequency  vector  for  three  angles  relative  to 
an  axis  of  the  detector.  For  the  blur-free  detector,  the  transfer 
function  is  simply  the  product  of  the  sines  in  the  two  direc¬ 
tions  induced  by  the  binning  operation.  The  OTF  of  the  blur- 
free  detector  is  nonzero  well  beyond  the  highest  frequency 
supported  by  the  lattice.  Any  component  of  an  input  signal  at 
these  higher  frequencies  will  contribute  to  a  lower  frequency 
alias  in  the  output,  as  per  Eq.  (38),  and  while  it  is  not  obvi¬ 
ous  from  the  point  of  view  of  frequency  space  the  sum  over 
aliases  in  Eq.  (38)  will  be  precisely  equivalent  to  the  detector 


Wiener  Spectra 


Frequency  Ip/mm 


Fig.  3.  The  Wiener  spectra  lV(f)/|/l|  for  the  three  model  detectors  as  in  Fig. 
2.  The  residual  additive  white  noise  <J>|A|/m  +  We/|A|  has  been  set  to  0. 


Fig.  4.  GDQE  as  a  function  of  frequency  for  the  three  model  detectors  as  in 
Fig.  2,  with  the  residual  white  noise  4>|A|/m  +  W£/|A|  set  to  0. 


simply  binning  each  incident  x-ray.  For  the  alias-free  detec¬ 
tors,  there  is  very  little  response  to  frequencies  beyond  those 
supported  by  the  lattice.  For  each  detector,  three  angles  are 
plotted,  but  the  angular  dependence  for  the  alias-free  detec¬ 
tors  is  small  enough  to  not  be  apparent  on  the  graph. 

The  Wiener  spectra  are  shown  in  Fig.  3.  Again,  instead  of 
plotting  a  function  of  two  variables,  fx  and  fy ,  we  plot  the 
NPS  as  a  function  of  the  magnitude  of  the  frequency  vector 
for  three  angles,  0=0,  27,  and  45°  (27°  corresponds  to  a 
slope  of  1:2  relative  to  the  lattice).  However,  the  NPS  shows 
little  angular  dependence.  For  the  Wiener  spectrum  one  only 
needs  to  look  at  frequency  values  supported  by  the  lattice, 
i.e.,  /,e[-  1/2L,1/2L]  and  fy e [ -  1/2L,1/2L]  (for  conve¬ 
nience  one  can  consider  W  to  be  periodic  in  the  frequency 
plane).  Thus,  at  6=0°  one  only  needs  to  graph  up  to  1/2Z, 
=  3.5  mm-1,  but  at  0=45°  the  frequencies  are  on  the  diago¬ 
nal  of  the  square,  so  one  goes  up  to  \2I2I .  =  4.9  min  *.  At 
0=27°,  one  goes  up  to  1/(2L  cos  0)=4mirTl.  For  this  graph 
the  constant  offset  d>|A|/w  +  W£-/|A|  has  been  set  to  zero. 
For  the  blur-free  detector,  the  NPS  is  flat,  which  follows 
mathematically  from  Eq.  (89)  and  the  fact  that  /(  f)  for  these 
detectors  is  simply  related  to  sine  functions,  or  more  physi¬ 
cally  by  noting  that  for  a  detector  which  simply  bins  incident 
x-rays  adjacent  cells  will  be  uncorrelated  so  the  NPS  is  flat. 
For  the  alias-free  detectors,  the  NPS  is  suppressed  by  factors 
of  the  square  of  the  transfer  function. 

The  GDQE  as  a  function  of  frequency  are  shown  in 
Figs.  4-6  for  a  range  of  values  of  the  residual  white  noise 
$>\A\lm  +  WEl\A\.  For  each  graph,  values  are  plotted  as  a 
function  of  the  magnitude  of  the  spatial  frequency  for  angles 
0,  27,  and  45°  relative  to  an  axis.  In  each  case,  the  GDQE 
falls  off  most  quickly  at  0=0°  and  least  quickly  at  0=45°, 
which  represents  the  fact  that  on  the  diagonal  the  sampling 
rate  is  increased  by  a  factor  of  vX.  In  the  case  where  the 
residual  white  noise  is  zero  (Fig.  4),  the  GDQE  of  the  blur- 
free  detector  drops  like  the  square  of  a  sine  function.  For  the 
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GDQE 


SNR2  as  a  function  of  displacement 


(50|i  wire  20  mm  long  parallel  to  axis) 


Displacement  (mm) 


Fig.  5.  GDQE  as  a  function  of  frequency  for  the  three  model  detectors  as  in 
Fig.  2,  with  the  residual  white  noise  &\A\lm  +  WEl\A\  set  to  42. 


alias-free  detectors,  the  GDQE  remains  at  nearly  unity  up 
to  the  lattice  cutoff,  the  factor  of  T2(  f)  canceling  the  same 
factor  in  the  colored  part  of  the  noise.  The  GDQE  of  the 
blur-free  detector  shows  some  response  beyond  the  lattice 
cutoff.  Though  small,  this  portion  of  the  GDQE  is  physical 
and  it  will  be  shown  that  the  responses  to  the  aliased 
frequencies  can  not  be  trivially  dismissed.  With  the  addition 
of  residual  white  noise  the  GDQE  of  all  three  models  is 
reduced,  as  shown  in  Fig.  5  ($>\A\lm+WEl\A\  =  42)  and 
Fig.  6  (<I>|A|/m  +  WeI\A\~$2).  These  figures  illustrate  that 
the  alias-free  detectors  are  more  sensitive  to  sources  of  re¬ 
sidual  white  noise  than  blur-free  detectors.  Indeed,  in  Fig.  6 
the  blur-free  detector  now  has  higher  GDQE  than  the 
\= 0.463  mm  detector  even  at  low  spatial  frequencies.  From 
the  spatial  point  of  view  this  is  quite  reasonable.  A  detector 


GDQE 


<D|A|/m  +  WE/|Aj  =  82 


Fig.  6.  GDQE  as  a  function  of  frequency  for  the  three  model  detectors  as  in 
Fig.  2,  with  the  residual  white  noise  4>|A|/m  +  We!\A\  set  to  82. 


Fig.  7.  Relative  SNR2  as  a  function  of  position  for  a  50  wide,  20  mm 
long  wirelike  object  parallel  to  one  axis  of  the  detector.  The  horizontal  axis 
of  the  graph  gives  the  displacement  of  the  wire,  so  that  at  0  mm  the  wire  is 
over  a  single  column  of  sensitive  elements,  while  at  0.07  mm  the  wire 
straddles  two  columns.  The  verticle  scale  is  arbitrary  (dependent  upon  the 
contrast  of  the  wire). 


whose  transfer  function  is  designed  to  remove  aliases  has  a 
relatively  wide  point  spread  function,  and  therefore  a  rela¬ 
tively  wide  autocovariance  function.  The  ideal  observer 
makes  use  of  digital  values  in  array  elements  whose  distance 
from  the  position  of  the  signal  is  up  to  several  times  the 
lengths  of  the  tails  of  these  functions,  so  for  the  same  input 
signal  the  ideal  observer  will  have  to  integrate  over  a  larger 
region  on  an  alias-free  detector  and  thus  will  be  more  sensi¬ 
tive  to  any  residual  uncolored  noise. 

While  GDQE  is  directly  related  to  the  average  value  of 
SNR2  by  Eq.  (82),  high-frequency  signals  (i.e.,  x-ray  pro¬ 
jection  images  of  small  objects  or  objects  whose  projected 
density  varies  quickly  with  position)  can  demonstrate  signifi¬ 
cant  changes  in  SNR2  with  position.  To  explore  this,  con¬ 
sider  the  SKE/BKE  task  associated  with  an  object  50  ^m 
wide  and  20  mm  long.  Wires  of  this  width  have  been  used  in 
neurological  and  cardiovascular  stents.49  Figure  7  shows  the 
SNR2  for  such  an  object,  oriented  parallel  to  an  axis  of  the 
detector,  as  a  function  of  displacement  in  the  direction  of  the 
shorter  (50  /nm )  axis.  The  scale  of  the  vertical  axis  is  arbi¬ 
trary  as  we  won’t  set  the  inherent  contrast  of  the  signal.  In 
Fig.  7,  the  0  mm  displacement  corresponds  to  the  signal 
being  centered  over  the  sensitive  elements  of  the  detectors. 
For  the  blur-free  detector,  at  0  mm  displacement  the  signal 
falls  into  a  single  column  of  detectors,  so  the  SNR2  corre¬ 
sponds  to  counting  statistics  for  one  column  of  elements  20 
mm  long.  The  SNR2  is  constant  until  the  0.05  mm  mark, 
after  which  the  signal  is  shared  between  two  columns  of 
detectors,  resulting  in  a  drop  in  SNR2.  Physically,  the  num¬ 
ber  of  x-rays  attenuated  by  the  object  is  independent  of  its 
position,  but  at  a  displacement  of  0.7  mm  the  signal  is  shared 
equally  between  two  columns,  and  since  the  noise  is  as¬ 
sumed  to  be  uncorrelated,  the  variance  in  the  total  counts  in 
the  two  columns  is  twice  that  of  one  column,  so  the  SNR2  is 
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Table  I.  SNR2  averaged  over  position  and  orientation  for  the  projection 
image  of  an  object  0.05  by  20  mm.  As  the  SNR2  scales  with  the  square 
of  the  contrast  of  the  image,  only  relative  values  are  meaningful.  The  ± 
represent  the  rms  (root-mean-square)  fluctuations  in  the  SNR2  with  position, 
not  the  statistical  uncertainties. 


SNR2  for  50  /*m  wide,  20  mm  long  object 

<t>\A\lm  +  WEl\A\=0 

=  42 

=  82 

Blur-free 

1.81  ±0.02 

1.79  ±0.02 

1.73  ±0.02 

Alias-free  (X— 0.463  mm) 

2.16±0.16 

1.87  ±0.07 

1.48  ±0.02 

Alias-free  (X=0.34  mm) 

2.15±0.15 

1.99±0.10 

1.68  ±0.05 

reduced  by  half.  The  alias-free  detectors  show  less  sensitivity 
to  position,  as  the  signal  is  always  shared  between  multiple 
columns.  As  before,  curves  are  shown  for  three  values  of  the 
residual  white  noise  <J)\A\/m  +  WEi\A\(02A2i  and  8~).  For 
all  models,  the  SNR2  drops  as  the  residual  white  noise  in¬ 
creases,  but  this  effect  is  greater  for  the  alias-free  models. 
Table  I  gives  the  average  SNR2  for  detection  of  the  0.05 
wide,  20  mm  long  wire,  now  averaged  over  both  position 
and  orientation.  Additionally,  the  root-mean-square  variation 
in  SNR2  is  given,  to  indicate  the  degree  to  which  the  detect¬ 
ability  of  the  wire  would  vary.  Again,  the  alias-free  detectors 
give  a  higher  SNR2  if  the  residual  white  noise  is  kept  suffi¬ 
ciently  low.  In  calculating  the  SNR2  of  the  projection  of  the 
wire  using  Eq.  (78),  if  the  summation  over  aliases  is 
dropped,  the  resulting  integral  decreases  by  about  10%.  Thus 
the  contributions  of  the  aliased  signal  to  the  SNR  are  not 
always  negligible. 

Somewhat  speculatively  one  can  consider  tasks  which  de¬ 
pend  upon  higher  frequency  components  of  the  signal/ 
Consider  a  5  mm  square  with  10%  contrast,  and  a  second 
square  whose  edges  have  been  smoothed  by  convolving  with 
a  0.15  mm  rect  function,  so  that  the  resulting  signal  “ramps 
up”  over  a  distance  of  0.3  mm.  The  SNR2  for  the  SKE/BKE 
task  of  distinguishing  between  these  two  objects  is  given  in 
Table  II.  It  is  interesting  to  note  that,  mathematically,  the 
SNR2  is  sufficiently  large  that  the  ideal  observer  can  perform 
this  task  efficiently,  although  whether  a  human  could  do  this 
is  questionable.  On  the  other  hand,  edge  detection  is  impor¬ 
tant  both  computationally  and  probably  as  part  of  the  strat¬ 
egy  of  human  observers,  so  the  ability  to  perform  this  task  is 
not  a  priori  irrelevant.  Again,  the  alias-free  detectors  have  a 
higher  SNR2  if  the  residual  white  noise  is  zero,  but  as  the 
task  now  depends  more  heavily  on  the  higher  frequency 

Table  II.  SNR2  averaged  over  position  and  orientation  for  distinguishing 
between  the  projection  of  a  5  mm  square  with  sharp  boundaries  and  a  square 
whose  boundaries  ramp  up  over  a  0.3  mm  region.  Normalization  corre¬ 
sponds  to  a  10%  contrast  with  an  x-ray  flux  corresponding  to  1400  x-rays/ 
pixel.  The  ±  represents  the  rms  fluctuations  in  the  SNR2  with  position,  not 
the  statistical  uncertainties. 


SNR2  for  discontinuity  vs  slope 


<h|A|/m  +  W£/|A|=0 

=  42 

=  82 

Blur-free 

14±2 

14±2 

13±2 

Alias-free  (X =0.463  mm) 

16±3 

12±  1 

7.6±0.5 

Alias-free  (X=0.34  mm) 

16±3 

n 
+  1 

10±  1 
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Table  III.  SNR2  averaged  over  position  and  orientation  for  detecting  the 
projection  of  a  0.5  mm  square.  Normalization  is  arbitrary.  The  ±  represents 
the  rms  fluctuations  in  the  SNR2  with  position,  not  the  statistical  uncertain¬ 
ties. 


SNR2  for  0.5  mm  square 

Blur-free 

Alias- free 

4>|A|/m  +  VMA|=0 

1.91  ±0.01 

2.08  ±0.01 

=  42 

1.89±0.01 

1.99±0.01 

=  82 

1.83±0.01 

1.789±0.003 

(X =0.463  mm) 
Alias-free 

2.07±0.01 

2.01  ±0.01 

1.855  ±0.006 

(X=0.34  mm) 

components  of  the  signal,  the  alias-free  detectors  are  more 
sensitive  to  residual  white  noise,  with  the  crossover  at 
<$>\A\lm  +  WE/\A\=42.  The  detection  ofa0.5  mm  square  is 
shown  in  Table  III.  Here,  the  lower  (but  nonzero)  frequen¬ 
cies  dominate  the  response  of  the  detector,  so  that  in  general 
the  antialiasing  detectors  gain  from  the  removal  of  the 
aliased  noise  without  losing  any  signal. 

It  is  interesting  to  note  that  the  SNRs  for  the  tasks  and 
models  described  above  do  not  vary  greatly.  Many  factors 
not  considered  here  will  greatly  effect  the  performance  of 
real  detectors,  beginning  with  the  fact  that  less  than  100%  of 
the  incident  x-rays  will  produce  secondary  quanta.  The  color 
of  the  Wiener  spectrum  need  not  be  the  same  as  the  transfer 
function,  due  to,  for  example,  x-rays  interacting  at  various 
depths  in  the  detector.51  The  efficiency  of  collection  of  the 
secondary  quanta  can  also  have  significant  effects.52  For  Csl 
detectors,  of  which  our  “alias-free”  detector  is  a  rough 
model,  the  fill  factor  is  a  minor  effect  so  long  as  the  ampli¬ 
fication  m  is  sufficiently  large.48  For  selenium  detectors,  of 
which  our  “blur-free”  detector  is  an  approximation,  it  is 
possible  to  have  an  effective  fill  factor  significantly  greater 
than  the  geometric  fill  factor  of  the  TFT  array.53  54  In  any 
case,  our  purpose  here  is  merely  to  indicate  some  of  the 
issues  which  must  be  faced  in  quantifying  digital  detectors  of 
these  types.  In  addition,  we  did  not  consider  geometric  fac¬ 
tors  such  as  x-ray  focal  spot  size  and  x-ray  parallax48  which 
reduce  the  high-frequency  content  of  the  incoming  signals. 

VIII.  DISCUSSION  AND  CONCLUSION 

The  results  of  this  paper  set  up  a  framework  for  quantita¬ 
tive  measurements  of  digital  systems  in  a  manner  analogous 
to  the  now  common  analysis  of  screen-film  systems  in  terms 
of  the  gray-scale  transfer  characteristic,  the  optical  transfer 
function,  the  Wiener  spectrum,  and  signal-to-noise  ratio.  The 
logic  of  this  framework  is  by  design  close  to  the  classic  work 
on  film-screen  systems.  While  many  pieces  of  this  argument 
have  appeared  in  the  work  of  others,  as  noted  throughout  the 
text,  it  seemed  desirable  to  produce  a  coherent  systematic 
exposition.  The  results  can  be  seen  as  appropriate  limits  of 
those  of  Barrett  et  al„  but  the  theoretical  construction  here 
emphasizes  the  parallels  with  the  classic  results  on  screen- 
film  systems.  While  detectors  consisting  of  discrete  elements 
do  not  have  continuous  translational  symmetry,  the  remain¬ 
ing  discrete  symmetry  allows  one  to  use  the  appropriate  Fou¬ 
rier  technique.  The  advantage  of  this,  as  in  the  screen-film 
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case,  is  that  one  can  explicitly  solve  for  the  mask  function  of 
the  ideal  observer  [see  Eq.  (75)]  and  thereby  obtain  the  op¬ 
timal  SNR.  As  for  screen-film,  this  formula  can  be  inter¬ 
preted  as  the  ratio  of  the  square  of  the  signal  in  each  fre¬ 
quency  bin  to  the  noise  in  each  frequency  bin,  as  measured 
by  the  Wiener  spectrum,  integrated  over  bins. 

We  have  investigated  several  highly  idealized,  but  not 
completely  unrealistic,  models  of  detectors,  and  illustrated 
some  of  the  issues  inherent  in  various  design  decisions.  This 
analysis  is  incomplete  and  intended  to  point  toward  issues 
which  could  be  addressed  in  other  work.  However,  our  re¬ 
sults  suggest  that  for  typical  tasks  the  detectability  of  objects 
as  determined  by  SNR2  is  not  drastically  affected  by  the 
decision,  in  and  of  itself,  to  suppress  or  not  to  suppress 
aliases.  In  any  real  device,  of  course,  this  design  decision  is 
linked  to  many  other  parameters.  This  article  should  be  of 
use  in  clarifying  what  is  actually  experimentally  measured  in 
testing  such  devices. 

The  results  presented  here  are  exact  for  the  SKE/BKE 
task  as  approached  by  the  ideal  observer  under  the  assump¬ 
tions  of  linearity,  homogeneity,  and  stationarity.  However, 
each  of  these  assumptions  is  only  approximately  true  in  prac¬ 
tice.  The  finite  extent  of  real  detectors  trivially  shows  that 
they  are  not  homogeneous,  but  for  a  variety  of  tasks  edge 
effects  are  negligible.  More  importantly,  many  digital  detec¬ 
tors  in  practice  show  significant  inhomogeneity  and  nonsta- 
tionarity.  The  work  of  Barrett  et  al  is  sufficiently  general  to 
cover  these  cases.  Further,  the  inhomogeneity  and  nonsta- 
tionarity  of  a  given  instrument  often  occur  in  ways  which  are 
different  for  each  individual  device,  so  that  while  the  extra 
information  is  relevant  to  the  particular  device  one  has  mea¬ 
sured,  the  extra  information  is  often  not  generalizable  to 
other  devices  of  the  same  manufacture.  This  extra  informa¬ 
tion  is  useful  for  optimizing  certain  tasks  using  the  particular 
device,  but  of  less  use  in  understanding  a  class  of  devices. 
Additionally,  while  the  signal  detection  task  of  the  ideal  ob¬ 
server  under  SKE/BKE  conditions  certainly  shares  some  fea¬ 
tures  with  the  task  which  human  observers  face,  and  has 
under  many  conditions  been  shown  to  correlate  well  with  the 
ability  of  human  observers  (for  example,  screen-film  images 
of  nylon  beads55),  it  is  still  a  very  idealized  task.  For  ex¬ 
ample,  if  edge-detection  is  important,  then  higher  frequency 
parts  of  the  incoming  signal  become  more  important  than 
would  be  expected  given  simply  the  signal  detection  task. 
From  the  point  of  view  of  a  radiologist,  a  clear  edge  might  be 
used  to  identify  and  distinguish  the  existence  of  a  lesion 
from  a  variation  in  projected  density  of  the  underlying  organ, 
particularly  in  the  presence  of  the  “structured  noise”  of 
other  anatomical  features. 

Clearly  these  are  issues  for  further  study,  but  while  OTF, 
W,  SNR,  and  GDQE  are  certainly  useful  because  they  are 
objectively  measurable  and  in  a  mathematically  precise  man¬ 
ner  are  related  to  tasks  which  approximate  those  of  the  hu¬ 
man  observer,  it  is  worth  remembering  that  measurement  of 
these  quantities  does  not  obviate  the  need  for  observer  stud¬ 
ies,  particularly  with  practicing  radiologists. 
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A  method  is  proposed  for  generating  synthetic  mammograms  based  upon  simulations  of  breast 
tissue  and  the  mammographic  imaging  process.  A  computer  breast  model  has  been  designed  with  a 
realistic  distribution  of  large  and  medium  scale  tissue  structures.  Parameters  controlling  the  size  and 
placement  of  simulated  structures  (adipose  compartments  and  ducts)  provide  a  method  for  consis¬ 
tently  modeling  images  of  the  same  simulated  breast  with  modified  position  or  acquisition  param¬ 
eters.  The  mammographic  imaging  process  is  simulated  using  a  compression  model  and  a  model  of 
the  x-ray  image  acquisition  process.  The  compression  model  estimates  breast  deformation  using 
tissue  elasticity  parameters  found  in  the  literature  and  clinical  force  values.  The  synthetic  mammo¬ 
grams  were  generated  by  a  mammogram  acquisition  model  using  a  monoenergetic  parallel  beam 
approximation  applied  to  the  synthetically  compressed  breast  phantom.  ©  2002  American  Asso¬ 
ciation  of  Physicists  in  Medicine.  [DOI:  10. 1 1 1 8/1 . 1 50 1 143] 
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image  acquisition 


I.  INTRODUCTION 

Visibility  of  breast  lesions  in  mammography  is  compromised 
by  overlapping  projections  of  normal  anatomic  structures 
that  generate  a  background  texture,  which  can  mask  existing 
abnormalities  or  introduce  false  ones.  Several  authors1"3 
have  shown  that  these  parenchymal  patterns  are  often  the 
limiting  factor  in  detection  tasks.  A  3D  simulation  of  mam¬ 
mography  is  proposed  to  provide  insight  into  the  formation 
of  such  patterns.  The  simulation  allows  one  to  analyze  the 
correlation  between  the  3D  composition  of  the  breast  and  its 
2D  mammographic  appearance.  By  identifying  the  dominant 
anatomical  structures  found  in  an  average  breast,  this  model 
can  help  analyze  the  deformation  of  those  structures  during 
the  exam  and  their  appearance  in  mammograms.  A  3D  mam¬ 
mography  simulation  can  also  serve  as  a  complement  to  ex¬ 
periments  with  respect  to  positioning,  compression,  and  ac¬ 
quisition.  Optimization  of  imaging  parameters  (such  as 
compression  angle  and  force,  x-ray  tube  kVp  and  mAs,  etc.) 
cannot  be  achieved  by  repeatedly  imaging  the  same  patient, 
due  to  concerns  about  the  radiation  dose.  Such  simulations 
can  also  be  useful  in  training  medical  personnel  by  demon¬ 
strating  the  effects  of  technique  selection  on  image  quality  or 
by  determining  3D  lesion  position  from  two  or  more  projec¬ 
tions.  Finally,  a  3D  breast  model  can  provide  a  theoretical 
framework  for  testing  new  breast  imaging  modalities.  Many 
new  modalities  are  being  developed  today,  including 
stereoscopy,4-6  tomosynthesis.7  and  3D  image  recon¬ 
struction,8  which  are  expected  to  provide  more  diagnostic 
information  about  normal  and  abnormal  tissue  structure.  It  is 


essential  in  the  development  of  such  systems  to  have  a  tool 
that  can  be  used  to  test  the  visibility  of  breast  structures  and 
help  select  optimal  views  for  3D  reconstruction,  since  the 
number  of  views  is  limited  by  the  amount  of  radiation  re¬ 
ceived. 

Historically,  mammography  simulation  started  with  the 
design  of  the  first  mathematical  breast  models  for  computing 
the  dose  received  by  a  patient  during  an  examination  using 
Monte  Carlo  simulation  of  x-ray  interactions.9  These  simu¬ 
lations  have  used  fairly  crude  models  of  breast  anatomy, 
lacking  internal  structures.  More  recent  analytical  models  of 
mammographic  image  acquisition  have  related  the  average 
values  of  the  incident  x-ray  flux,  linear  attenuation  coeffi¬ 
cients  of  breast  tissue,  and  the  film  density  or  pixel  digital 
values  in  the  obtained  mammograms.10 

There  are  two  approaches  to  modeling  the  image  content 
of  mammograms.  In  a  2D  approach,  mammograms  are  mod¬ 
eled  based  upon  the  analysis  of  spatial  correlation  between 
image  pixel  values,  using  various  random  field  methods.11-13 
Such  models  can  match  some  of  the  statistical  properties  of 
real  mammograms,  but  they  cannot  reveal  the  relationship 
between  the  3D  structures  of  the  breast,  nor  they  can  consis¬ 
tently  produce  images  of  the  same  breast  with  modified  po¬ 
sition  or  acquisition  parameters. 

In  this  paper  we  propose  a  second  approach,  whereby 
mammograms  are  modeled  by  projection  of  simulated  3D 
anatomic  structures,  based  upon  the  size  and  the  distribution 
of  large  and  medium  scale  tissue  regions  found  in  the  breast. 
It  is  our  hypothesis  that  the  distribution  of  the  2D  structures 
seen  in  mammograms  reflects  the  distribution  of  the  3D  tis- 
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Breast  Tissue  Model 


Large  Scale  Elements:  AT,  FGT 
Medium  Scale  Elements:  Shells,  Blobs , 
andCDuctalTree) 


Synthetic  Mammograms 


Fig.  1.  Components  of  a  system  for  mammography  simulation.  The  flow¬ 
chart  shows  the  order  (top-to-bottom)  in  which  the  model  components  are 
simulated. 

sue  structures  of  the  breast.  Thus,  the  background  texture  in 
the  synthetic  mammograms  so  generated,  should  have  simi¬ 
lar  properties  to  those  found  in  clinical  images. 

Positioning  and  compression  significantly  affect  the  ap¬ 
pearance  of  mammograms.  Until  recently14”16  there  were  no 
models  of  breast  deformation,  due  to  the  complex  anatomy 
of  the  different  types  of  interwoven  breast  tissue,  whose  me¬ 
chanical  properties  are  difficult  to  analyze.  We  have  approxi¬ 
mated  breast  compression  by  separate  deformations  of  tissue 
layers  positioned  normal  to  the  compression  plates. 

The  objective  of  this  work  is  to  generate  synthetic  mam¬ 
mograms.  A  method  for  achieving  this  goal  is  described  in 
Sec.  II,  while  the  results  of  the  simulation  are  shown  in  Sec. 
III.  An  accompanying  paper  details  an  analysis  of  the  quality 
of  the  model.17 

II.  MAMMOGRAPHY  SIMULATION 

The  proposed  mammography  simulation  consists  of  three 
major  components:  a  3D  software  breast  phantom,  a  com¬ 
pression  model,  and  an  x-ray  image  acquisition  model  (see 
Fig.  1).  The  breast  phantom  is  a  software  tissue  model  con¬ 
taining  two  ellipsoidal  regions  of  large  scale  tissue  elements: 
predominantly  adipose  tissue  (AT)  and  predominantly  fibro- 
glandular  tissue  (FGT).  The  internal  tissue  structures  of  these 
regions,  namely  the  adipose  compartments  and  the  breast 
ductal  network,  are  approximated  by  realistically  distributed 
medium  scale  phantom  elements:  shells,  blobs,  and  the  simu¬ 
lated  ductal  tree.  The  compression  model  is  based  upon  tis¬ 
sue  elasticity  properties  and  a  breast  deformation  model.  De¬ 
formation  is  simulated  separately  for  tissue  slices  (or  layers) 


(a) 


(b) 

Fig.  2.  Anatomic  structures  of  the  breast  included  in  the  tissue  model,  (a) 
La ige  scale  regions  seen  in  MRI  images  as  predominantly  adipose  tissue 
region.  AT  (bright),  and  predominantly  fibroglandular  region,  FGT  (dark 
surrounded  by  the  AT),  (b)  Subgross  (thick)  histologic  slice  showing  large 
scale  regions:  AT  and  FGT,  and  medium  scale  tissue  structures:  compart¬ 
ments  surrounded  by  Cooper’s  ligaments  in  the  AT  and  small  adipose  cavi¬ 
ties  and  ducts  within  the  FGT.  (Subgross  breast  histology  image  provided 
courtesy  of  Dr.  R.  D.  Cardiff.) 

positioned  normal  to  the  compression  plates.  Each  slice  is 
approximated  by  a  beam  composed  of  two  different  tissues. 
Deformed  slices  are  stacked  to  produce  a  model  of  the  com¬ 
pressed  breast.  The  mammogram  acquisition  model  was 
adopted  from  the  literature,10  assuming  monoenergetic  x  rays 
and  a  parallel  beam  geometry  without  scatter.  Presently,  the 
synthetic  mammograms  are  generated  with  a  spatial  resolu¬ 
tion  of  200  /xm/pixel  because  the  current  model  version  does 
not  include  fine  tissue  details.  This  resolution  is  comparable 
to  the  resolution  of  digitized  mammograms  in  the  Mini 
MIAS  database  (obtained  by  averaging  4X4  pixels  in  the 
original  MIAS  database18). 

A.  Software  breast  phantom 
1.  Modeling  large  scale  tissue  regions 

Figure  2  illustrates  the  types  of  anatomic  structures  of  the 
breast  which  are  included  in  the  breast  phantom.  Figure  2(a) 
is  an  MRI  breast  section,  showing  the  shape  and  position  of 
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Fig.  3.  Orthogonal  sections  of  the  uncompressed  breast  model.  Parameters 
without  subscripts  correspond  to  the  semi-axes  of  the  ellipsoidal  approxima¬ 
tion  of  the  breast  outline;  superscripts  A  and  P  correspond  to  the  anterior  and 
posterior  border  of  the  FGT  model  region,  respectively. 


(a)  (b)  (c) 


Fig.  4.  Examples  of  different  tissue  distributions  in  real  breasts,  illustrated 
with  clinical  mammograms  from  the  MIAS  database.  The  amount  of  adipose 
tissue  affects  the  size  and  visibility  of  compartments  seen  in  the  AT  and  FGT 
regions. 


The  interiors  of  the  shells  and  blobs  have  the  elastic  and 
x-ray  attenuation  properties  of  adipose  tissue.  As  a  first  ap¬ 
proximation,  the  adipose  compartments  are  represented  by 
spheres.  The  size  of  the  spheres  can  vary  to  allow  for  normal 
breast  anatomic  variations. 

Figure  3  shows  two  orthogonal  cross  sections  of  a  breast 
tissue  phantom.  Simulated  regions  of  predominantly  adipose 
and  predominantly  fibroglandular  tissue  are  seen,  together 
with  the  spherical  approximation  of  adipose  compartments  in 
those  regions.  Note  that  the  size  of  the  simulated  compart¬ 
ments  in  the  AT  and  FGT  regions  differ.  The  size  of  the 
adipose  compartments  varies  in  different  women,  depending 
upon  the  amount  of  adipose  tissue  in  the  breast,  as  seen  in 


the  large  scale  tissue  regions.  The  MRI  shows  the  predomi¬ 
nantly  fibroglandular  tissue,  appearing  as  a  dark  region  in  the 
central  part  of  the  breast,  and  the  predominantly  adipose  tis¬ 
sue,  a  brighter  region  surrounding  the  fibroglandular  tissue. 
Internal  structures  of  these  tissue  regions  are  less  clearly  vis¬ 
ible,  due  to  the  relatively  low  MRI  resolution  of  approxi¬ 
mately  1  mm/pixel.  Figure  2(b)  shows  a  subgross  histologic 
slice  of  the  breast,  obtained  after  mastectomy.  The  predomi¬ 
nantly  fibroglandular  region  in  the  slice  is  represented  by  a 
darker  image  region  in  the  center,  surrounded  by  brighter, 
predominantly  adipose  tissue.  Adipose  tissue  is  organized 
into  round  compartments,  formed  by  fibrous  Cooper's  liga¬ 
ments.  The  FGT  region  also  contains  adipose  compartments, 
but  they  are  smaller  in  size  than  the  compartments  in  the  AT 
region.  Analysis  of  subgross  histologic  slices  of  the  breast 
and  the  corresponding  mammograms  showed  that  projec¬ 
tions  of  these  compartments  dominantly  contribute  to  the 
formation  of  parenchymal  patterns.  Therefore,  simulated  adi¬ 
pose  compartments  were  included  as  medium  scale  breast 
model  elements. 

2.  Modeling  adipose  tissue  compartments 

The  adipose  compartments  are  approximated  by  thin 
shells  in  the  AT  region  and  small  blobs  in  the  FGT  region. 


(a)  (b)  (c) 


Fig.  5.  Sections  of  breast  models  with  different  sized  tissue  elements  labeled 
according  to  Table  I  as:  (a)  “Small,”  (b)  “Medium,”  and  (c)  “Large.” 
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Fig.  6.  Examples  of  computer  generated  ductal  lobes,  (a)  A  simulated  mam¬ 
mogram  with  five  duct  lobes.  For  the  purpose  of  this  illustration,  other 
medium  scale  image  elements  have  been  suppressed,  (b)  Two  views  of  the 
same  simulated  ductal  network,  used  to  generate  image  in  (a).  (The  letter 
“NT  indicates  position  of  the  nipple.) 


Fig.  4.  Adipose  compartments  are  more  easily  identified  in  a 
histologic  slice  than  in  a  mammogram,  since  the  latter  image 
contains  the  superimposed  projections  of  many  tissue  layers. 
Simulations  of  breasts  with  different  sized  adipose  compart¬ 
ments  are  illustrated  in  Fig.  5.  The  procedure  for  generating 
simulated  tissue  compartments  starts  by  filling  the  3D  breast 
phantom  with  compartments  of  the  largest  selected  size  until 
they  start  to  intersect  each  other.  The  compartment  size  is 
then  reduced  and  the  procedure  continued  until  the  smallest 
selected  size  has  been  reached. 

3.  Modeling  breast  ductal  network 

To  achieve  a  sufficiently  realistic  tissue  phantom  it  is  also 
necessary  to  model  the  breast  ductal  network.  Ducts  can  be 
visualized  by  galactography,  a  clinical  x-ray  imaging  proce¬ 
dure  whereby  the  ducts  are  enhanced  by  injection  of  a  con¬ 
trast  agent.  The  breast  ductal  system  consists  of  about  15-20 
ductal  lobes,  each  corresponding  to  a  major  duct  branching 
from  the  nipple  into  a  network  of  smaller  ducts.  In  galacto- 
grams,  usually  a  single  lobe  is  enhanced.  Larger  ducts  are 
more  visible  than  the  smaller  ones,  since  they  attenuate  more 
x  rays. 

The  main  focus  of  our  model  is  on  the  pattern  of  duct 
branching.  This  pattern  can  be  expressed  by  a  ramification 
matrix,  representing  probabilities  of  branching  at  different 
levels  of  a  tree  structure.19  There  are  no  previous  reports  in 
the  literature  on  analyzing  the  breast  ductal  network  by  rami¬ 
fication  matrices.  The  class  of  random  binary  trees  was  cho¬ 


sen  for  modeling  the  ducts  since  it  offers  the  least  con¬ 
strained  branching  pattern.19  The  ductal  model  consists  of  15 
lobes  and  each  lobe  is  simulated  by  a  different  random  bi¬ 
nary  tree.  Different  trees  are  generated  by  using  different 
random  number  generator  seeds.  A  simulated  ductal  network 
is  shown  in  Fig.  6.  For  clarity,  only  5  out  of  15  lobes  are 
visible. 


B.  Mammographic  compression  model 

The  mammographic  compression  simulation  is  based 
upon  a  deformation  model  including  realistic  tissue  elasticity 
properties.  Elasticity  parameters  of  the  tissue  found  in  the 
literature  vary  significantly.20-22  One  of  the  reasons  for  these 
variations  is  that  the  experiments  for  determining  the  elastic 
properties  have  been  performed  using  small  samples  taken 
from  a  particular  tissue  type  (e.g.,  adipose,  fibroglandular,  or 
cancerous).  However,  the  breast  is  comprised  of  a  compli¬ 
cated  admixture  of  different  tissues  which  affect  the  elastic 
behavior  of  the  whole  organ. 

The  parameters  most  often  used  for  description  of  elastic 
properties  are  the  Young’s  linear  elasticity  modulus,  E ,  and 
the  Poisson’s  ratio,  v,  defined  by 


<t  FIA  _  Aw 

£=7=a/7t  v~~K. T' 


(i) 


The  Young’s  modulus  relates  the  strain,  e,  as  a  measure  of 
deformation  (i.e.,  the  fractional  change  in  length  A///)  and 
the  stress,  o*  (the  force  F  applied  to  the  surface  area  A  of  the 
deformed  object).  The  Poisson  ratio,  v ,  is  equal  to  the  ratio 
of  transverse  contraction,  Aw,  to  the  elongation.  A/,  of  a 
deformed  bar.  It  is  usually  assumed  that  human  tissue  can  be 
approximated  as  an  incompressible  material,  whose  volume 
does  not  change  during  deformation.20  For  incompressible 
materials  v^0.5. 

The  other  elasticity  moduli  used  to  describe  the  behavior 
of  material  are  the  bulk  modulus,  K ,  and  shear  modulus,  G , 
which  can  be  expressed  using  the  values  of  E  and  v\ 


K= 


3(l-2v)9 


G~ 


2(1  +  v) 


(2) 


There  is  a  relationship  between  the  bulk  elasticity  modulus, 
K ,  material  density,  p,  and  the  speed  of  sound  through  the 
material,  u ,  given  by 

v~\[K/p.  (3) 

In  the  compression  model,  the  elasticity  parameter  values  of 
adipose  and  fibroglandular  tissues  were  computed  using  the 
values  of  ultrasound  velocity  through  various  tissues  found 
in  the  literature.23 

The  Mammography  Quality  Standards  Act24  regulates  the 
minimum  and  maximum  breast  compression  to  be  used  in 
mammography.  In  general,  the  mammography  technician 
will  apply  the  maximum  force  tolerated  by  the  patient  to 
achieve  optimum  quality  mammograms.  Sullivan  et  aL25  re¬ 
ported  a  statistical  analysis  of  the  compression  force  and 
compressed  breast  thickness,  measured  during  560  exams. 
Simulated  force  values  were  selected  to  address  these  con¬ 
siderations. 
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Fig.  7.  Separate  deformation  of  breast  model  slices.  Il¬ 
lustrated  are  the  steps  of  rectangular  slice  approxima¬ 
tion,  deformation  of  the  approximated  slice,  and  gen¬ 
eration  of  the  final  shape  of  the  deformed  slice.  These 
steps  are  followed  by  stacking  all  deformed  slices  to¬ 
gether  to  generate  the  compressed  breast  shape  (not 
shown). 


There  are  very  few  published  results  about  deformation  of 
breast  tissue  structures  during  mammography.  One  reason  is 
that  it  is  hard  to  provide  3D  visualization  of  such  deforma¬ 
tions.  Clinically  available  3D  breast  imaging  techniques,  ul¬ 
trasound  and  MRI,  have  lower  resolution  and  poorer  quality 
than  mammography.  In  addition,  both  methods  use  much  less 
compression  and  are  not  applicable  for  analyzing  tissue  de¬ 
formation  during  mammography.  In  our  compression  model, 
tissue  deformation  is  estimated  in  two  phases.  First,  the  large 
scale  model  elements,  the  AT  and  FGT  regions,  are  deformed 
to  determine  the  shape  of  the  compressed  breast.  Second,  the 
medium  scale  model  elements  are  deformed  by  transforming 
the  shells  and  spheres  into  ellipsoids.  After  compression  and 
x-ray  image  acquisition,  the  medium  scale  elements  appear 
as  elliptical  structures  in  synthetic  mammograms  corre¬ 
sponding  to  the  oval  shaped  lucencies  seen  in  real  mammo¬ 
grams. 

The  compression  of  large  scale  model  elements  is  simu¬ 
lated  in  the  form  of  separate  deformations  of  the  breast  tissue 
slices,  positioned  normal  to  the  compression  plates.  In  this 
paper,  the  medio-lateral  oblique  (MLO)  mammographic  view 
is  modeled  since  it  provides  visualization  of  more  breast  tis¬ 
sue  than  other  views.26  Moreover,  in  a  number  of  countries, 
e.g.,  UK,  the  Netherlands,  and  Sweden,  it  is  the  only  view 
obtained  by  breast  screening.27  Compression  for  other  views, 
e.g.,  cranio-caudal  (CC),  can  be  simulated  by  a  modification 
of  the  described  model.  In  the  case  of  MLO  compression, 
tissue  slices  are  positioned  normal  to  the  MLO  view  plane. 
In  this  case,  a  slice  appears  as  a  semiellipse,  and  the  FGT 
portion  of  the  slice  appears  as  the  intersection  of  two  semiel¬ 
lipses,  as  shown  in  Fig.  3.  Slice  thickness  corresponds  to  the 
model  resolution  specified  at  the  beginning  of  the  simulation. 


Deformation  of  each  slice  is  computed  in  three  steps.  First,  a 
rectangular  slice  approximation  is  computed.  Second,  slice 
deformations  are  estimated  using  a  composite  beam  model. 
Third,  the  compressed  slice  shape  is  computed  from  the  de¬ 
formed  rectangular  approximation.  Compressed  slices  are 
stacked  together  to  form  the  compressed  breast  model  shape. 
A  detailed  description  of  the  processing  steps  is  given  in  the 
following. 

1.  Rectangular  slice  approximation 

A  slice  of  the  breast  model  is  replaced  by  its  rectangular 
approximation.  The  whole  slice  region  and  its  FGT  portion 
are  approximated  by  rectangles,  satisfying  the  following  con¬ 
straints:  (i)  the  area  of  the  rectangular  slice  approximation, 
A  snce>  and  area  of  the  rectangular  FGT  approximation,  AFG, 
are  the  same  as  the  corresponding  areas  in  the  original  slice; 
(ii)  the  side  of  the  rectangular  slice  approximation  in  the 
chest-nipple  direction,  </SUcc ,  and  the  side  of  the  rectangular 
FGT  approximation  in  the  same  direction,  dFG ,  are  equal  to 
the  corresponding  dimensions  in  the  original  slice;  and  (iii) 
the  distance  between  the  centers  of  gravity  of  the  rectangular 
slice  approximation  and  the  rectangular  FGT  approximation 
is  the  same  as  the  corresponding  distance  in  the  original 
slice.  This  constraint  is  not  used  for  deformation  of  the  slice 
through  the  nipple,  because  it  may  produce  a  rectangular 
approximation  of  the  FGT  region  protruding  outside  of  the 
slice.  Instead,  the  FGT  region  is  positioned  so  that  it  touches 
the  nipple  side  of  the  rectangular  approximation  of  the  whole 
slice. 

Using  these  constraints,  the  dimensions  and  relative  posi¬ 
tions  of  the  rectangular  slice  and  FGT  approximation  are 
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computed.  An  example  of  a  slice  through  the  nipple  is  illus¬ 
trated  in  Fig.  7. 


Ea  ^  PA 

E?g  *fg  pFGUpG 


(7) 


2.  Slice  deformation  via  a  composite  beam  model 

The  approximating  rectangles  preserve  the  elastic  proper¬ 
ties  of  the  corresponding  slice  regions.  The  elastic  properties 
of  the  AT  and  FGT  regions  are  modeled  by  linear  Young’s 
moduli,  Ea  and  EFG ,  respectively.  This  rectangular  slice  ap¬ 
proximation  can  be  treated  as  a  composite  2D  elastic  beam, 
positioned  between  two  bars  corresponding  to  slices  through 
the  compression  plates.  Deformation  is  estimated  by  apply¬ 
ing  a  force  to  the  compression  plates  which  in  turn  deforms 
the  composite  beam.  It  is  assumed  that  (i)  the  slice  thickness 
is  much  smaller  than  the  sides  of  the  approximating  rect¬ 
angles  and  (ii)  during  the  compression,  the  slice  stays  in  the 
same  plane  as  before  compression.  The  latter  assumption  is 
based  on  the  fact  that  in  reality  slices  of  tissue  are  not  com¬ 
pressed  independently;  the  neighboring  tissue  (above  or  be¬ 
low)  partially  confine  the  slice  to  a  plane.  Slices  in  the  planes 
above  or  below  the  nipple  level  do  not  have  initial  contact 
with  the  compression  plates.  It  is  assumed  that  all  slices  de¬ 
form  with  the  same  strain  value  which  is  equal  to  the  strain 
of  the  slice  in  the  nipple  level: 


*Slicc=4iPcf= 


A  W  RecSlice 
W  RecSlice 


=  1- 


wncw 
w  RecSlice 

^RecSlice 


1 


^Compress 
w  Relax  ’ 

(4) 


where  €SIicc  and  e^Ic  represent  strain  for  any  slice  and 
strain  for  the  slice  in  the  nipple  level,  respectively;  w^slice 
and  w RecSlice  represent  size  of  the  rectangular  approximation 
normal  to  the  compression  plates,  before  and  after  compres¬ 
sion,  respectively;  wRelax  and  vvCompress  represent  breast  thick¬ 
ness  before  and  after  compression,  respectively.  The  inten¬ 
sity  of  the  compression  force  is  included  indirectly,  by 
specifying  the  thickness  of  the  compressed  breast.  Since  a 
linear  model  of  tissue  elasticity  is  used,  only  the  ratio  of 
Young’s  moduli  for  the  AT  and  FGT  region  is  needed  to 
estimate  slice  deformation. 

The  stress  in  the  rectangular  FGT  region,  aR^  is  the  same 
as  the  stresses  crR  and  crR  in  the  parts  of  the  rectangular  AT 
region  with  thickness  wR[  and  wR 3, 


(Tr]  ”  °rR2~(rRy 


(5) 


Replacing  cry  by  ,  from  Eq.  (1),  the  strain  in  the  rectan¬ 
gular  FGT  approximation  is  calculated  as: 


Ea 


“  ^FG  1 


(6) 


The  ratio  EA  fEFG  is  computed  based  upon  the  relationship 
between  the  bulk  elastic  moduli,  tissue  density,  and  velocity 
of  sound  propagation  through  the  tissue.  Measured  values  of 
the  velocity  of  sound  in  samples  of  adipose  and  fibroglandu- 
lar  tissue  are  va=1470  m/s  and  i;FG=  1545  m/s, 
respectively.23  Densities  of  the  adipose  and  fibroglandular 
tissue  are  pA  =  930  kg/m3  and  pFG=  1040  kg/m3, 
respectively.28  Using  Eq.  (3),  gives 


Finally,  there  is  a  relationship  between  the  strains  in  dif¬ 
ferent  parts  of  the  rectangular  slice  approximation,  since: 

AwL=Awi?i  +  Awy?2+Awi?3,  (8) 

which  after  dividing  by  wL  and  using  the  fact  that  €R^  —  €R^ 
due  to  symmetry,  yields 


—  €r 


1  WL 


+  €j? 


Wr2 

1  WL  * 


(9) 


Equations  (4),  (6),  and  (9)  yield  dimensions  of  the  rectan¬ 
gular  slice  and  FGT  region  normal  to  the  compression  plates 
after  the  compression,  w^slice  and  w^c¥G,  respectively: 

<ZSnce=”d  1  -  «l),  wr^fG=w*2(  1  -  %)•  (10) 


Assuming  that  the  areas  of  the  rectangular  slice  and  FGT 
region,  ASUce  and  AFG,  stay  the  same  before  and  after  com¬ 
pression,  we  can  compute  the  dimension  of  the  rectangular 
slice  and  FGT  approximation  in  the  chest-nipple  direction 
after  compression,  d^slice  and  dJf™FG ,  respectively: 


new  _  ^  Slice  mew  ^FG 

^RecSlice  new  ’  ^RecFG  new 

w  RecSlice  WRecFG 


(ID 


3 .  Compressed  slice  from  the  deformed 
rectangular  approximation 

The  final  step  of  the  slice  deformation  modeling  is  the 
computation  of  the  compressed  breast  slice  from  its  de¬ 
formed  rectangular  approximation.  The  compressed  breast 
does  not  have  an  ellipsoidal  but  rather  a  flattened  shape.“6'"9 
The  thickness  of  the  compressed  breast  is  constant  and  equal 
to  the  distance  between  the  compression  plates,  vvComprcss, 
everywhere  except  in  a  narrow  region  close  to  the  front  edge 
of  the  breast.  Analysis  of  that  region  on  a  mammogram  was 
used  to  estimate  the  breast  thickness  directly  from 
mammograms.30  To  achieve  a  realistic  shape  of  the  com¬ 
pressed  breast  slice,  a  correction  was  applied  to  the  model.31 
This  correction  assumes  that  the  deformed  breast  slice  con¬ 
sists  of  a  rectangle  positioned  at  the  chest  wall  side,  and  a 
semiellipse  attached  to  the  rectangle,  extending  forward  to 
the  nipple  (see  Fig.  7). 

Parameters  of  the  deformed  rectangle  and  semiellipse  are 
computed  satisfying  the  following  constraints:  (i)  the  sum  of 
the  rectangular  area,  ARectanglc,  and  the  semielliptical  area, 
A  Ellipse^  js  equa]  to  the  area  0f  the  whole  uncompressed  slice, 
Aslice;  (ii)  one  side  of  the  rectangle  and  one  axis  of  the 
semiellipse  are  equal  to  the  distance  between  the  compres¬ 
sion  plates  for  the  compressed  breast,  wComPrcss;  and  (iii)  the 
slice  region  where  the  thickness  is  less  then  vvComprcss  con¬ 
tains  10%  of  the  whole  mammogram  breast  area.30  The  de¬ 
scribed  correction  for  flattening  the  compressed  breast  is 
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(a)  (b)  (c) 


Fig.  8.  Examples  of  synthetic  images  with  different  sizes  of  tissue  model 
elements  (i.e.,  simulated  compartments  in  the  AT  and  FGT),  which  are  used 
for  comparison  with  clinical  mammograms.  The  images  are  labeled  accord¬ 
ing  to  Table  I  as:  (a)  “Small/’  (b)  “Medium,”  and  (c)  “Large.” 


only  used  to  determine  the  border  of  the  whole  compressed 
slice.  Deformation  of  the  FGT  region  is  still  computed  using 
the  2D  composite  beam  approximation. 

Separate  processing  of  individual  model  slices  is  followed 
by  stacking  the  deformed  slices  together  to  get  the  3D  com¬ 
pressed  breast  model.  Slices  whose  relaxed  (noncompressed) 
thickness  is  less  than  the  compressed  breast  thickness  are  not 
processed  at  all;  they  are  assumed  to  preserve  their  relaxed 
shape.  The  compressed  breast  thickness  was  used  instead  of 
the  compression  force  to  compute  the  deformation  of  the 
breast  model  slices.  The  compression  force  can  be  calculated 
from  the  difference  of  the  relaxed  and  compressed  breast 
thickness  and  tissue  elastic  moduli  using  Hooke’s  law.  When 
computing  the  force,  the  values  for  both  the  elastic  moduli  of 
the  adipose  and  fibroglandular  tissue  are  needed. 

C.  X-ray  mammogram  acquisition  model 

The  x-ray  image  acquisition  model  consists  of  an  x-ray 
propagation  model,  which  includes  attenuation  by  the  breast 
tissue  and  conversion  of  the  x-ray  energy  into  film  density, 
and  a  model  of  mammographic  film  digitization.  In  the  case 
of  digital  mammography,  the  film  and  digitization  models 
should  be  replaced  by  a  model  of  a  solid-state  x-ray  detector 
array.  The  model  is  adopted  from  the  literature,10  and  for 
simplicity  assumes  a  monoenergetic  x-ray  spectrum  and  a 
parallel  beam  geometry,  without  scatter. 

The  mammogram  acquisition  model  relates  the  spatial 
distribution  of  the  x-ray  energy  imparted  to  the  intensifying 
screen,  £/,  to  the  mammogram  digital  values  (after  digitiza¬ 
tion),  DV,  and  the  linear  x-ray  attenuation  coefficient  of  tis¬ 
sue,  fit : 


{f*  -  =  vv  Compress 

-J  fi,(x,y,z)dz },  (12) 


DV(x,y)=a-by\ogl0{(3EI(x,y)}, 


(13) 


Table  I.  Radii  of  simulated  adipose  compartments  in  the  AT  and  FGT,  used 
to  generate  the  synthetic  mammograms. 


Structure  radii 

Model  regions 

AT  (mm) 

FGT  (ram) 

Small 

2.7-6.7 

1.3-2.7 

Medium 

4-10 

2-4 

Large 

5.3-13.3 

2.7-S.3 

where  rj  represents  the  quantum  efficiency  of  the  screen,  <f>  is 
the  fluence  at  the  entrance  to  the  breast,  E  the  x-ray  photon 
energy,  C  is  the  attenuation  factor  due  to  the  compression 
paddle  and  grid,  wCompress  is  the  compressed  breast  thickness, 
a  and  b  are  digitization  coefficients,  and  y  and  ft  are  the  film 
gamma  and  the  speed  coefficient,  respectively.  It  is  assumed 
that  the  digitization  output  is  proportional  to  optical  density, 
thus  Eqs.  (12)  and  (13)  can  be  simplified  so  that  DV  is  lin¬ 
early  proportional  to  the  ray  sum.  The  linear  x-ray  attenua¬ 
tion  coefficients  of  the  AT  and  FGT  tissue,  taken  from  the 
literature,32  are  /xAr=  0.456  cm-1  and  /xFG=  0.802  cm-1,  at 
20  keV.  Further  details  of  the  acquisition  model  are  given  in 
the  literature.33 

III.  SIMULATION  RESULTS  AND  DISCUSSION 

Figure  8  shows  three  synthetic  mammograms  generated 
by  our  simulation.  These  three  images  differ  in  the  size  of  the 
simulated  medium  scale  tissue  structures.  The  ranges  of  adi¬ 
pose  compartments  radii  simulated  in  the  mammograms 
(shown  in  Fig.  8)  are  given  in  Table  I. 

It  can  be  noted  that  the  proportions  of  the  breast  model, 
i.e.,  vertical  to  horizontal  dimension  ratio  after  compression, 
agree  with  the  “standard  breast”  from  Novak,29  defined  by 
averaging  dimensions  of  27  compressed  breasts.  Dimensions 
of  the  breast  model  are  smaller  than  the  “standard  breast”  by 
approximately  15%. 

The  synthetic  mammograms  were  printed  on  film  (AGFA 
LR5200,  AGFA-Gevaert,  Belgium)  lifesize,  and  were  shown 
to  radiologists  in  the  Breast  Imaging  Center  of  Thomas  Jef¬ 
ferson  University.  Qualitatively,  subregions  of  synthetic  and 
clinical  images  were  reported  to  have  similar  appearance 
when  viewed  at  a  distance  of  1-2  m.  When  examined 
closely,  it  was  observed  that  the  synthetic  images  lack  blood 
vessels  and  other  organized  fine  tissue  structures.  In  addition, 
the  borders  between  the  AT  and  FGT  regions  in  the  synthetic 
mammograms  appeared  as  a  clear,  geometrically  regular 
separation  degrading  the  subjective  perception  of  reality.  For 
the  latter  reason,  the  model  was  modified  by  addition  of 
small  random  variations  to  the  position  of  the  borders  be¬ 
tween  the  compressed  AT  and  FGT  phantom  regions.  This 
correction  is  included  in  Fig.  8. 

As  stated  in  Sec.  I,  the  model  for  producing  the  synthetic 
mammograms  was  based  upon  the  hypotheses  that  the  size 
and  the  distribution  of  simulated  3D  tissue  elements  are  simi¬ 
lar  to  those  found  in  the  real  breasts,  and  that  the  3D  tissue 
distribution  is  reflected  in  the  distribution  of  2D  mammo¬ 
graphic  structures.  In  order  to  evaluate  the  synthetic  images, 
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we  performed  statistical  comparisons  of  several  texture  de¬ 
scriptors  computed  in  synthetic  and  clinical  mammograms, 
including:  average  size  of  image  objects,  texture  energy,  and 
fractal  dimension.  These  descriptors  have  been  previously 
used  in  the  literature  for  the  analysis  of  parenchymal  patterns 
in  mammograms.34-36  The  image  objects’  size,  analyzed  us¬ 
ing  mathematical  morphology,  can  be  related  to  the  size  of 
the  simulated  3D  tissue  structures.  Texture  energy  and  fractal 
dimension  are  sensitive  to  small  scale  changes  in  image  in¬ 
tensities,  corresponding  to  fine  tissue  detail.  Details  about  the 
texture  analysis  of  synthetic  and  clinical  mammograms  are 
given  in  the  accompanying  paper.17  Quantitatively,  the  syn¬ 
thetic  mammograms  have  a  similar  distribution  of  values  av¬ 
eraged  over  a  large  number  of  clinical  mammograms.  The 
best  matching  was  observed  for  the  synthetic  images  gener¬ 
ated  using  the  simulated  adipose  compartments  with  radii  of 
4-13.3  mm  in  the  AT  regions,  and  radii  of  2.7-5.33  mm  and 
1. 3-2.7  mm  in  the  retroareolar  and  dense  FGT  regions, 
respectively.17  It  is  expected  that  the  introduction  of  detailed 
tissue  structures  in  the  breast  model  will  enhance  the  local 
variations  of  synthetic  mammograms  and  the  variations  in 
feature  distribution  needed  to  better  match  real  images.  In 
addition,  the  simulated  ducts  and  the  compression  model 
were  separately  evaluated  and  compared  with  clinically  ac¬ 
quired  data.33,37 

IV.  CONCLUSIONS 

A  method  is  described  for  generating  synthetic  mammo¬ 
grams  using  simulations  of  breast  tissue  and  the  mammo- 
graphic  imaging  process.  A  software  breast  phantom  was  de¬ 
veloped,  which  contains  realistic  large  and  medium  scale 
tissue  structures,  derived  from  an  understanding  of  the  mac¬ 
roscopic  anatomic  tissue  organization.  Parameters  control¬ 
ling  the  size  and  placement  of  the  tissue  simulating  structures 
provide  flexibility  to  generate  a  large  database  of  synthetic 
images  with  different  characteristics.  Mammographic  imag¬ 
ing  is  simulated  using  a  compression  model  and  a  model  of 
the  x-ray  image  acquisition.  The  compression  model  esti¬ 
mates  breast  deformation  using  tissue  elasticity  parameters 
found  in  the  literature  and  realistic  values  of  compression 
force.  The  synthetic  mammograms  were  generated  by  a 
mammogram  acquisition  model  using  a  monoenergetic  par¬ 
allel  beam  approximation,  applied  to  the  synthetically  com¬ 
pressed  breast  phantom. 

The  proposed  simulation  can  be  used  in  analysis  of  breast 
positioning,  compression,  and  image  acquisition  parameters. 
The  software  breast  phantom  can  be  used  as  a  test  object  for 
optimizing  mammographic  systems  or  testing  novel  systems 
for  3D  reconstruction  of  breast  images.  The  simulation  can 
be  used  for  analyzing  the  correlation  between  the  3D  breast 
composition  and  2D  mammogram  characteristics,  e.g.,  the 
parenchymal  patterns,  which  can  be  used  for  estimation  of 
the  cancer  risk.38  Computer  algorithms  for  characterization 
of  normal  breast  tissue  could  be  tested  on  large  databases  of 
normal  images  with  random  variations  of  tissue  structures, 
generated  by  the  model.  Synthetic  mammograms  with  simu¬ 
lated  abnormalities  or  abnormalities  extracted  from  clinical 


mammograms  could  be  used  for  testing  algorithms  for  can¬ 
cer  detection. 
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We  have  evaluated  a  method  for  synthesizing  mammograms  by  comparing  the  texture  of  clinical 
and  synthetic  mammograms.  The  synthesis  algorithm  is  based  upon  simulations  of  breast  tissue  and 
the  mammographic  imaging  process.  Mammogram  texture  was  synthesized  by  projections  of  simu¬ 
lated  adipose  tissue  compartments.  It  was  hypothesized  that  the  synthetic  and  clinical  texture  have 
similar  properties,  assuming  that  the  mammogram  texture  reflects  the  3D  tissue  distribution.  The 
size  of  the  projected  compartments  was  computed  by  mathematical  morphology.  The  texture  energy 
and  fractal  dimension  were  also  computed  and  analyzed  in  terms  of  the  distribution  of  texture 
features  within  four  different  tissue  regions  in  clinical  and  synthetic  mammograms.  Comparison  of 
the  cumulative  distributions  of  the  mean  features  computed  from  95  mammograms  showed  that  the 
synthetic  images  simulate  the  mean  features  of  the  texture  of  clinical  mammograms.  Correlation  of 
clinical  and  synthetic  texture  feature  histograms,  averaged  over  all  images,  showed  that  the  syn¬ 
thetic  images  can  simulate  the  range  of  features  seen  over  a  large  group  of  mammograms.  The  best 
agreement  with  clinical  texture  was  achieved  for  simulated  compartments  with  radii  of  4-13.3  mm 
in  predominantly  adipose  tissue  regions,  and  radii  of  2.7-5.33  and  1.3 -2.7  mm  in  retroareolar  and 
dense  fibroglandular  tissue  regions,  respectively.  ©  2002  American  Association  of  Physicists  in 
Medicine.  [DOI:  10.1118/1.1501144] 

Key  words:  mammography  simulation,  3D,  synthetic  mammograms,  texture  analysis 


I.  INTRODUCTION 

We  have  proposed  an  approach  to  generate  synthetic  mam¬ 
mograms  based  upon  a  3D  simulation  of  mammography.1 
Synthetic  mammographic  texture  is  produced  by  projecting 
simulated  3D  breast  anatomic  structures.  In  clinical  images, 
the  overlapped  projections  of  normal  anatomic  tissue  struc¬ 
tures  generate  a  background  texture  in  mammograms  which 
can  mask  the  existing  abnormalities  or  introduce  false  ones. 
The  simulation  can  be  used  to  optimize  positioning,  com¬ 
pression  and  acquisition  in  order  to  improve  the  visibility  of 
the  breast  tissue,  and  to  test  new  breast  imaging  modalities. 

The  proposed  mammography  simulation  consists  of  three 
major  components.  First,  a  3D  software  breast  phantom  con¬ 
tains  two  ellipsoidal  regions  of  large  scale  tissue  elements: 
predominantly  adipose  tissue  (AT)  and  predominantly  fibro¬ 
glandular  tissue  (FGT)  regions.  Internal  structures  of  these 
regions,  namely  the  adipose  compartments  and  breast  ductal 
network,  are  approximated  by  realistically  distributed  me¬ 
dium  scale  phantom  elements:  shells  filled  with  simulated 
adipose  tissue  and  a  synthetic  ductal  tree.  Second,  a  com¬ 
pression  model  of  the  breast  deformation  occurring  during  a 
mammographic  exam  is  based  upon  tissue  elasticity  proper¬ 
ties.  Deformation  is  simulated  separately  for  layers  of  tissue 
positioned  normal  to  the  compression  plates.  Each  slice  is 


approximated  by  a  rectangular  beam  composed  of  AT  and 
FGT  regions.  The  slices  are  computationally  deformed,  as¬ 
suming  clinical  values  of  the  compression  force.  Deformed 
slices  are  stacked  together  to  produce  a  model  of  the  com¬ 
pressed  breast.  Third,  mammogram  image  acquisition  is 
modeled  assuming  monoenergetic  x  rays  and  a  parallel  beam 
geometry  without  scatter.  Details  of  the  simulation  are  given 
in  the  accompanying  paper.1 

Ideally,  each  of  the  three  components  of  the  simulation 
should  be  evaluated  separately  by  a  3D  imaging  technique. 
There  is,  however,  a  significant  difference  in  tissue  properties 
captured  by  the  clinically  available  3D  breast  imaging  mo¬ 
dalities  (ultrasound  and  MRI)  and  mammography  which  is 
the  focus  of  our  simulation.  Breast  ultrasound  and  MRI  also 
have  different  resolution  and  compression  than  mammogra¬ 
phy.  With  these  issues  in  mind  we  have  evaluated  the  tissue 
model  indirectly,  assuming  that  a  relationship  exists  between 
the  distribution  of  3D  breast  tissue  structures  and  the  2D 
parenchymal  pattern.  It  is  our  hypothesis  that  the  texture 
properties  computed  in  synthetic  and  clinical  images  have 
similar  distributions. 

There  are  two  approaches  to  mammogram  synthesis 
found  in  the  literature:  (i)  direct  modeling  of  2D  distribution 
of  pixels  and  (ii)  simulation  of  3D  tissue  distribution  and  the 
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mammographic  imaging.  Bochud  et  al?  modeled  mammo¬ 
gram  texture  as  a  “clustered  lumpy  background”  by  random 
placement  of  “blob”  clusters,  visually  resembling  tissue  ap¬ 
pearance  in  mammograms.  Synthetic  images  were  evaluated 
by  comparing  their  power  spectra  and  statistical  moments 
with  the  values  from  32  clinical  mammograms.  Good  agree¬ 
ment  of  the  first  and  the  second  moments  in  clinical  and 
synthetic  images  were  observed,  with  similar  statistical  prop¬ 
erties  overall.  Heine  et  al?  modeled  a  mammogram  as  evolv¬ 
ing  from  a  process  of  passing  a  random  field  (colored  noise) 
through  a  linear  filter  with  a  self-similar  characteristic,  based 
upon  the  analysis  of  60  clinical  mammograms.  Such  an  ap¬ 
proach  can  match  some  of  the  statistical  properties  of  clinical 
images  but  cannot  relate  the  3D  tissue  structures  and  their 
mammographic  appearance.  Both  papers  do  not  model  breast 
ducts  or  the  large  scale  tissue  regions.  Consequently,  the  im¬ 
ages  of  the  same  simulated  breast,  with  modified  positioning, 
compression,  or  x-ray  parameters  cannot  be  consistently  syn¬ 
thesized. 

Taylor  et  al?  generated  synthetic  images  by  mammogra¬ 
phy  simulation,  in  an  approach  similar  to  our  work.  The  fo¬ 
cus  of  their  simulation  is  on  modeling  breast  ducts  based 
upon  the  fractal  properties  of  the  duct  length  and  diameter. 
They  have  evaluated  the  synthetic  images  so  obtained  by 
comparing  the  Fourier  spectrum  with  that  computed  in  im¬ 
ages  of  tissue  slices  with  contrast  enhanced  ducts.  Good 
agreement  using  a  small  number  of  samples  was  observed. 

Separate  evaluations  were  performed  for  the  simulation  of 
the  ductal  network,  the  compression  model,  and  the  synthetic 
parenchymal  pattern.  Initial  feasibility  tests  of  the  ductal 
model  and  compression  simulation  are  presented 
elsewhere.5,6  This  paper  describes  the  analysis  of  the  syn¬ 
thetic  mammogram  texture. 

Synthetic  images  were  generated  by  simulating  the  x-ray 
image  acquisition  on  a  computationally  compressed  phan¬ 
tom.  Images  of  the  phantoms  were  generated  containing  dif¬ 
ferent  sizes  of  simulated  medium  scale  elements:  spherical 
shells  and  blobs.  The  synthetic  mammograms,  so  obtained, 
were  evaluated  by  comparing  them  with  clinical  images 
taken  from  the  MIAS  database  of  digitized  mammograms.7 
Subimages  taken  from  regions  corresponding  to  different  tis¬ 
sues  were  compared  separately,  including  the  subcutaneous 
AT,  retromammary  AT,  retroareolar  FGT,  and  dense  FGT  re¬ 
gions.  Three  texture  features  were  used  for  description  of  the 
parenchymal  pattern:  (i)  the  average  size  of  image  structures 


(a)  (b)  (c) 

Fig.  1.  Illustration  of  morphological  closing,  (a)  The  original  image  with 
objects  of  different  size,  (b)  The  image  of  the  structuring  element,  (c)  The 
resulting  image  obtained  by  the  morphological  closing  with  the  structuring 
element  from  (b)  applied  to  the  image  from  (a). 


computed  using  mathematical  morphology,  (ii)  the  texture 
energy,  and  (iii)  the  fractal  dimension.  Feature  values  were 
computed  over  each  clinical  and  synthetic  subimage  and  sta¬ 
tistically  compared  using  the  Kolmogorov-Smimov  test  and 
histogram  correlation.  Details  of  the  analysis  of  synthetic  and 
clinical  mammographic  texture  are  given  in  Sec.  II  and  the 
results  of  the  comparison  are  discussed  in  Sec.  III. 

II.  TEXTURE  ANALYSIS  OF  SYNTHETIC 
MAMMOGRAMS 

A.  Texture  descriptors 

The  following  texture  descriptors  were  used  for  the  evalu¬ 
ation  of  synthetic  mammogram  texture.  First,  size  analysis 
was  performed  by  a  sequence  of  morphological  closings  with 
disks  of  increasing  size  as  structuring  elements.8  Average 
image  brightness  increases  after  the  closing  operation.  The 
change  in  brightness  as  a  function  of  the  disk  radius  is  re¬ 
lated  to  the  size  distribution  of  radiolucent  (adipose)  areas  in 
the  mammograms.  Second,  texture  energy  analysis  was  per¬ 
formed  by  convolving  each  image  with  a  small  mask.9  Treat¬ 
ing  gray  scale  image  intensity  as  the  height  of  a  3D  object, 
this  mask  is  sensitive  to  local  roughness  of  the  image  sur¬ 
face.  Third,  fractal  dimension  was  computed  by  the  blanket 
box  counting  method  of  self-similarity  analysis.10 


1.  Morphological  analysis  of  image  structure  size 

Morphological  image  analysis  is  based  upon  the  shape  of 
image  objects  and  is  used  to  simplify  image  data  while  pre¬ 
serving  shape  characteristics.  The  theory  of  mathematical 
morphology  is  discussed  in  the  books  of  Matheron11  and 
Serra.8  An  application  oriented  tutorial  of  morphological  im¬ 
age  processing  is  given  by  Haralick.12 

Morphological  operations  are  performed  on  a  set  of  image 
pixels  using  a  second  set  of  pixels  called  the  structuring  el¬ 
ement.  Definitions  of  the  basic  operations  are  given  in  the 
Appendix.  The  opening  operation  is  used  for  size  analysis  of 
bright  objects,  and  closing  for  the  analysis  of  dark  objects. 
This  analysis  is  sensitive  to  the  radiolucent  areas  of  the 
mammogram,  corresponding  to  the  adipose  tissue  which  ap¬ 
pears  darker  than  the  surrounding  tissue.  X  rays  are  less  at¬ 
tenuated  by  adipose  tissue,  producing  greater  film  density 
than  connective  tissue. 

The  gray  scale  closing  first  replaces  each  pixel  with  the 
maximum  from  its  neighborhood  defined  by  the  structuring 
element  (a  disk).  The  original  values  are  then  recovered  for 
all  of  the  pixels,  except  for  those  from  regions  which  are 
both  darker  than  their  surroundings  and  smaller  than  the 
structuring  element.  As  an  illustration,  Fig.  1(a)  shows  an 
image  with  several  objects  of  different  size.  After  the  closing 
operation  with  the  structuring  element  from  Fig.  1(b),  the 
resulting  image  is  given  in  Fig.  1(c).  It  can  be  seen  that  dark 
objects  smaller  than  the  structuring  element  have  been  elimi¬ 
nated;  the  resulting  image  is  thus  brighter  than  the  original. 
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This  is  the  basis  for  morphological  size  analysis,  whereby 
the  change  in  average  image  brightness  (i.e.,  the  total  pixel 
sum  after  the  closing)  is  used  to  describe  the  size  distribution 
of  the  image  objects.  The  derivative  of  the  brightness  as  a 
function  of  size  shows  the  contribution  of  the  objects  equal 
in  size  to  the  structuring  element. 

Morphological  size  analysis  of  mammograms  has  been 
reported  previously  in  the  literature.13,14  Behrens  and 
Dengler13  reported  examples  of  applying  morphological  size 
analysis  at  global,  regional,  and  local  image  levels;  analysis 
of  calcifications  was  presented  as  a  local  processing.  Miller 
and  Astley 14  used  morphological  size  analysis  to  segment  the 
FGT  region  from  mammograms.  They  used  the  opening  op¬ 
eration  which  is  dual  to  closing;  it  replaces  the  bright  regions 
smaller  than  the  structuring  element  by  their  dark  surround¬ 
ing  pixels.  The  overall  image  brightness  is,  thus,  reduced. 
However,  the  authors  did  not  analyze  the  relationship  be¬ 
tween  morphological  feature  values  and  the  physical  proper¬ 
ties  of  the  anatomic  structures.  Our  research  represents  a 
novel  application  of  morphological  image  analysis  as  a  result 
of  using  the  simulated  3D  tissue  structures  to  synthesize  pa¬ 
renchymal  patterns. 

2.  Texture  energy  analysis 

Texture  energy  features  are  the  statistical  estimates  of  the 
outputs  from  a  filter  bank  implemented  in  the  form  of  local 
linear  transformations.  They  were  introduced  with  the  goal 
of  achieving  texture  segmentation  and  description  at  each 
image  pixel,  corresponding  to  a  hypothetical  low  level  func¬ 
tion  of  the  human  visual  system.9  The  filter  bank  consists  of 
small  2D  convolution  masks  whose  coefficients  are  com¬ 
puted  as  the  product  of  ID  masks  with  different  numbers  of 
zero  crossings.  Contrast  invariance  of  the  filter  outputs  is 
achieved  by  the  normalization  with  the  output  of  the  filter 
sensitive  to  the  average  local  image  intensity.  The  absolute 
values  or  variances  of  the  convolved  images  are  used  for 
analysis.  A  generalization  of  this  approach  can  include  a 
larger  set  of  local  linear  transformations,  and  the  estimation 
of  higher  order  moments  of  the  output  channel  histograms.15 

In  mammogram  processing  applications14,16  texture  en¬ 
ergy  was  usually  computed  using  a  single  or  a  few  convolu¬ 
tion  masks.  The  mask  sensitive  to  image  “ripple"  was  found 
to  be  the  most  efficient  in  segmenting  potentially  abnormal 
regions  in  mammograms,  a  task  which  is  related  to  the  local 
roughness  of  the  image  surface.  Texture  energy  features  have 
also  been  used  in  mammogram  registration.1718 

The  mask  coefficients  are  given  in  the  Appendix.  A  5X5 
“ripple"  convolution  mask,  R5R5  [Eq.  (A2)],  was  used,  cap¬ 
turing  local  roughness  of  the  image  surface.  The  absolute 
values  of  the  convolved  data  were  averaged  on  a  15X15 
window  and  normalized  by  the  “level"  mask,  L5L5  [Eq. 
(A2)],  providing  contrast  invariance. 

3.  Fractal  analysis 

Fractal  dimension  describes  self-similarity  of  image  prop¬ 
erties  at  different  spatial  scales.  It  is  common  to  perform 
fractal  analysis  on  the  area  of  the  image  surface,  obtained  by 


considering  the  pixel  values  as  local  surface  heights.  This 
area  is  related  to  the  roughness  of  the  image  texture.  A  com¬ 
plete  definition  of  the  fractal  dimension  of  image  surface 
area  is  given  in  the  Appendix. 

There  are  numerous  reports  in  the  literature  on  fractal 
analysis  of  mammograms.19-22  Caldwell  et  a/.19  analyzed  the 
fractal  dimension  of  various  parenchymal  patterns  and  the 
difference  between  the  fractal  dimensions  computed  over  the 
whole  image  and  within  a  region  near  the  nipple.  A  feature 
space  defined  by  these  two  fractal  features  was  segmented 
and  a  relatively  good  agreement  with  the  original  Wolfe 
classification23  was  observed.  Mammographic  calcifications 
have  been  segmented  using  a  variety  of  methods  for  comput¬ 
ing  fractal  dimension,  including  box  counting,20  iterated 
function  systems,21  and  fractal  Brownian  motion.22 

We  computed  fractal  dimension  by  the  blanket  algo¬ 
rithm.10  This  method  has  been  used  previously  in  the  detec¬ 
tion  of  calcifications  in  mammograms.20  The  fractal  dimen¬ 
sion  is  computed  for  each  pixel  by  analyzing  the  local  image 
surface  around  the  pixel.  A  15X15  window  was  selected, 
centered  on  each  pixel.  This  corresponds  to  the  nonlinear 
averaging  window  size  used  in  the  texture  energy  method.  A 
log-log  plot  of  Alocal(€)  is  generated  for  the  local  surface 
around  each  pixel.  The  local  fractal  dimension  value  Dlocal  is 
computed  as  the  slope  through  three  points  on  the  log-log 
plot,  corresponding  to  the  scale  parameter  values  of  e—2 ,  3, 
and  4  pixels. 

B.  Image  selection 

The  following  criteria  were  used  for  selection  of  the  clini¬ 
cal  and  synthetic  mammograms  to  be  used  for  comparison. 
First,  the  clinical  images  had  to  represent  normal  breast  tis¬ 
sue.  Second,  the  glandularity  seen  in  the  mammograms 
should  approximately  represent  the  average  breast  glandular¬ 
ity  (not  too  dense  and  not  predominantly  adipose).  Third, 
spatial  resolution  of  the  clinical  and  synthetic  mammograms 
should  be  matched.  The  clinical  images  were  selected  from 
the  MIAS  database7  of  digitized  mammograms  and  the  syn¬ 
thetic  mammograms  were  generated  for  varying  properties  of 
the  medium  scale  elements,  i.e.,  different  sizes  of  simulated 
adipose  compartments  in  the  AT  and  FGT  regions.  In  addi¬ 
tion,  the  comparison  was  repeated  for  the  same  set  of  clinical 
and  synthetic  mammograms  at  a  reduced  resolution.  The  im¬ 
ages  with  reduced  resolution  were  generated  by  averaging 
2X2  blocks  of  pixels  from  the  original  mammograms. 

1.  Clinical  mammograms 

Sixty-five  mammograms  from  the  Mini-MIAS  database  of 
clinical  mammograms  were  used,  having  a  spatial  resolution 
of  200  /un/pixel.  The  Mini-MIAS  database  was  obtained  by 
averaging  4X4  pixel  blocks  in  the  original  MIAS  mammo¬ 
gram  database.7  This  resolution  is  sufficient  for  the  evalua¬ 
tion  of  our  synthetic  mammograms  since  presently  they  do 
not  include  fine,  small  scale  tissue  detail.  The  selected  im¬ 
ages  represent  normal  cases  in  the  MIAS  database  with  the 
background  tissue  classified  as  “fatty-glandular."  As  the 
sizes  of  adipose  compartments  differ  for  various  tissue  re¬ 
gions,  up  to  four  25  mmX25  mm  subimages  per  mammo- 
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Fig.  2.  Tissue  regions  used  in  texture  analysis,  illustrated  on  a  clinical  mam¬ 
mogram  from  the  MIAS  database:  (1)  subcutaneous  adipose  tissue,  (2)  ret¬ 
romammary  adipose  tissue,  (3)  retroareolar  fibroglandular  tissue  (immedi¬ 
ately  posterior  to  the  nipple),  and  (4)  dense  fibroglandular  tissue. 


gram  were  selected  manually,  giving  a  total  of  219  sub¬ 
images,  from  the  following  regions  (see  Fig.  2):  (1)  subcuta¬ 
neous  fat;  (2)  retromammary  fat;  (3)  retroareolar  glandular 
tissue,  immediately  posterior  to  the  nipple;  and  (4)  dense 
glandular  tissue.  If  the  extent  of  a  tissue  region  could  not  be 
unambiguously  determined,  or  if  it  was  too  small  for  a  sub¬ 
image  window,  the  corresponding  tissue  sample  was  ex¬ 
cluded  from  analysis. 

2.  Synthetic  mammograms 

Synthetic  images  were  generated  at  a  spatial  resolution  of 
200  fx m/pixel,  matching  that  of  the  database.  Four  subimages 
per  synthetic  mammogram  were  selected  from  different  re¬ 
gions  in  the  same  manner  as  for  the  clinical  images.  The 
positions  of  the  subimages  were  determined  from  the  known 
extent  of  the  large  scale  model  elements,  the  AT  and  FGT 
regions.  Model  parameters  controlling  the  distribution  of  me¬ 
dium  scale  tissue  structures,  modeled  by  shells  in  the  AT  and 
spheres  in  the  FGT  regions,  were  varied  to  match  the  statis¬ 
tical  properties  of  real  images.  Three  groups  of  synthetic 
mammograms  were  tested.  The  groups  consisted  of  ten  syn¬ 
thetic  mammograms  each,  generated  randomly  using  the 
same  range  of  size  of  simulated  adipose  tissue  compart- 
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ments.  The  ranges  of  compartment  sizes  differed  between  the 
groups  by  30%  (see  Table  I  and  Fig.  8  in  the  accompanying 
paper1). 


C.  Statistical  comparison 

Two  methods  were  used  for  statistical  comparison  of  the 
texture  features.  First,  feature  histograms  were  computed  for 
each  subimage.  Synthetic  histograms  were  then  averaged 
over  all  subimages  of  the  same  tissue  type  and  were  com¬ 
pared  with  similarly  computed  clinical  histograms.  The  cor¬ 
relation  between  the  corresponding  clinical  and  synthetic  av¬ 
eraged  histograms  was  used  to  measure  how  well  the 
synthetic  images  approximated  the  clinical  images.  Next, 
mean  feature  values  (i.e.,  the  histogram  first  moments)  were 
computed  for  each  subimage.  Distributions  of  these  means 
for  all  subimages  of  the  same  tissue  type  were  then  analyzed 
and  compared  with  the  distributions  of  means  of  the  clinical 
images,  using  the  Kolmogorov-Smimov  (KS)  test.24  The 
maximum  difference  between  the  cumulative  distribution 
functions  (CDFs)  of  the  clinical  and  synthetic  mean  feature 
values  was  used  as  another  measure  of  quality  of  mammo¬ 
gram  synthesis.  In  both  methods  the  average  texture  features 
were  compared,  thereby  testing  the  ability  of  the  simulation 
to  match  the  average  properties  of  a  large  set  of  clinical 
mammograms,  rather  than  simulating  an  image  of  a  particu¬ 
lar  breast. 


1.  Analysis  of  feature  histograms 

As  a  measure  of  similarity  between  the  real  and  synthetic 
feature  distributions,  the  correlation  between  the  feature  his¬ 
tograms  was  calculated  for  each  of  the  clinical  and  synthetic 
subimages,  and  averaged  over  all  subimages  of  the  same 
tissue  type.  In  the  case  of  size  analysis,  the  correlation  was 
computed  between  the  brightness  gradient  (as  a  function  of 
the  structuring  element  radius)  of  clinical  and  synthetic  im¬ 
ages.  In  the  following  text,  these  derivative  values  are  re¬ 
ferred  to  as  the  “average  histogram  of  the  size  analysis  fea¬ 
ture.'’ 

The  coefficient  of  correlation,  /?,  between  the  real,  hR , 
and  synthetic,  hs ,  histograms  averaged  over  all  subimages 
(in  a  given  category)  is  computed  as: 


W0M0 
V2,[M/)]22,[  MO]*’ 


where  the  summation  runs  over  histogram  bins  i. 


(1) 


2.  Kolmogorov- Smirnov  (KS)  test 

The  KS  test  compares  two  random  distributions  based 
upon  the  maximum  difference  between  their  CDFs.24  It  be¬ 
longs  to  a  group  of  nonparametric  methods  which  make  no 
assumptions  about  the  types  of  distributions  used.  The  maxi¬ 
mum  difference  between  two  CDFs,  £>,  is  a  measure  of  the 
discrepancy  between  the  two  sets  of  samples.  Kolmogorov 
showed  that  for  two  sets  of  samples  with  the  same  parent 
distribution,  the  CDF  of  D  is  given  asymptotically  by:24 


2144  Bakic  eta!.:  Mammogram  synthesis  using  a  3D  simulation.  II 


2144 


Fig.  3.  Texture  energy  histograms  of  the  FGT  from  clinical  (left)  and  synthetic  (right,  primed)  mammograms,  (a)  Sample  subimage,  (b)  Image  of  texture 
energy  values,  (c)  Texture  energy  histogram  (normalized  for  the  range  of  feature  values). 


CO 

lim  P(Dm,„^z)  =  1-22 

m ,  n  — * 


Xexp 
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m  n  I 
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where  m  and  n  are  the  numbers  of  samples  in  the  two  sets, 
Dm%„  is  the  maximum  CDF  difference  for  the  given  number 
of  samples,  and  P  is  the  probability  that  Dm  n  is  less  than  a 
given  value  z>  The  level  of  significance,  a,  is  defined  as 

P(Dmjl>da)  =  a,  (3) 


where  da  is  the  critical  value  of  Dm  n  corresponding  to  the 
significance  a .  Thus,  the  observed  discrepancy  between  a 
CDF  drawn  from  clinical  mammograms  and  a  CDF  drawn 
from  simulated  mammograms  can  be  quantified  in  terms  of 


the  significance,  a,  which  is  the  probability  that  a  greater 
discrepancy  than  observed  would  occur  due  to  chance  alone. 
The  relationship  between  a  and  da  for  various  sample  sizes 
is  tabulated  in  several  textbooks.24,25 

The  CDFs  of  statistics  for  each  subimage  from  the  clinical 
and  synthetic  mammograms  were  compared.  In  the  texture 
energy  analysis  and  the  fractal  analysis,  for  each  subimage 
the  appropriate  feature  value  was  averaged  over  all  of  the 
pixels  in  the  subimage,  and  this  average  was  used  as  a 
sample  for  the  KS  test.  In  the  morphological  analysis,  for 
each  subimage  the  first  moment  of  the  brightness  gradient 
was  used  as  a  sample  value. 

Both  the  KS  test  and  the  histogram  correlation  show  how 
well,  on  average,  the  synthetic  images  can  approximate  the 
properties  of  the  clinical  mammographic  texture.  The  differ- 
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Fig.  4.  Size  analysis  of  the  FGT  from  clinical  (left)  and  synthetic  (right,  primed)  mammograms,  (a)  Sample  subimage,  (b)  Result  of  closing  with  a  10  pixel 
disk  structuring  element,  (c)  Result  of  closing  with  a  40  pixel  disk,  (d)  Change  in  brightness  (sum  of  all  pixels)  before  and  after  closing,  (e)  Gradient  of  the 
brightness. 


ence  between  the  two  methods  is  that  the  KS  test  compares 
the  mean  feature  values  averaged  over  each  subimage,  while 
the  histogram  correlation  takes  into  account  the  range  of  fea¬ 
ture  values  computed  locally  at  each  pixel. 

3.  Illustration  of  the  analysis 

An  illustration  of  the  histogram  analysis  is  given  in  Fig.  3 
by  the  texture  energy  features  computed  on  subimages  of 
retroareolar  glandular  tissue.  The  histogram  of  a  clinical  FGT 
subimage  is  shown  on  the  left  and  of  a  synthetic  subimage  on 
the  right.  Histograms  averaged  over  all  clinical  and  over  all 


synthetic  subimages  are  shown  in  Fig.  5(c).  Figure  4  illus¬ 
trates  the  analysis  of  object  size  distribution  for  the  retroar¬ 
eolar  glandular  tissue.  The  left-hand  side  shows  the  results 
for  the  clinical  FGT  and  the  right-hand  side  for  the  synthetic 
FGT.  The  upper  graphs  show  the  average  brightness  (offset 
for  the  brightness  of  the  original  image)  after  each  morpho¬ 
logical  closing  as  a  function  of  structuring  element  (disk) 
size  in  pixels.  Note  that  the  output  images  get  brighter  with 
increasing  disk  size,  as  seen  in  the  examples  of  the  images 
obtained  for  the  disk  radii  of  10  and  40  pixels.  The  graph  in 
Fig.  4(e)  shows  the  gradient  of  the  features  graphed  in  Fig. 
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Fig.  5.  Comparison  between  clinical  and  synthetic  images  of  retroareolar  fibroglandular  tissue:  CDFs  of  the  (a)  mean  texture  energy  feature  and  (b)  mean 
fractal  analysis  feature:  average  histograms  of  the  (c)  mean  texture  energy  feature  and  (d)  mean  fractal  analysis  feature. 


4(d).  The  gradient  is  used  in  place  of  the  feature  histogram 
for  the  size  analysis.  (Running  averages  of  the  gradients  are 
shown  for  clarity.)  Mean  feature  values,  whose  distributions 
were  analyzed  by  the  KS  test,  are  also  indicated  on  the 
graphs  in  Fig.  4(e).  For  texture  energy  and  fractal  dimension, 
the  mean  feature  values  were  computed  as  the  first  moments 
of  the  feature  histograms  of  each  subimage.  CDFs  for  all 
clinical  and  three  groups  of  synthetic  subimage  means  are 
shown  in  Fig.  6(b). 

Figure  5  shows  the  CDF  of  the  mean  feature  values  and 
the  average  histograms  for  texture  energy  and  fractal  dimen¬ 
sion  computed  for  the  samples  of  retroareolar  glandular  tis¬ 
sue  regions  in  clinical  and  synthetic  mammograms.  Figure  6 


gives  the  cumulative  distributions  of  the  first  moments  of  the 
brightness  gradients  and  the  average  brightness  gradients  for 
the  size  analysis  of  the  clinical  and  synthetic  samples  of  the 
subcutaneous  adipose  tissue  (left)  and  the  retroareolar  glan¬ 
dular  tissue  (right). 

III.  RESULTS  AND  DISCUSSION 

Results  of  the  synthetic  texture  evaluation  are  presented  in 
the  form  of  graphs  of  the  histogram  correlation  coefficients 
and  maximum  CDF  differences  computed  for  several  texture 
features.  Simulations  were  performed  for  three  ranges  of 
compartment  size  and  for  two  spatial  resolutions.  Figures 
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Fig.  6.  Comparison  between  clinical  and  synthetic  images  using  size  analysis:  CDFs  of  the  mean  features  for  (a)  subcutaneous  adipose  tissue  and  (b) 
retroareolar  fibroglandular  tissue:  average  histograms  for  (c)  subcutaneous  adipose  tissue  and  (d)  retroareolar  fibroglandular  tissue. 


7-10  show  the  results  for  subcutaneous  adipose  tissue,  ret¬ 
romammary  adipose  tissue,  retroareolar  glandular  and  dense 
glandular  tissue  regions,  respectively.  The  abscissa  in  these 
graphs  is  not  a  continuous  variable,  but  indicates  the  size 
range  of  the  simulated  breast  anatomic  structures  (adipose 
compartments).  Texture  features  are  labeled  by  different 
symbols:  circle = average  structure  size,  diamond = texture 
energy,  and  triangle = fractal  dimension. 

The  graphs  labeled  (a)  in  Figs.  7-10  show  the  values  of 
histogram  correlation.  Higher  values  indicate  better  correla¬ 
tion,  with  a  maximum  possible  value  of  unity.  A  boot-strap 
analysis  of  the  histogram  correlation  (averaged  over  all  tissue 
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types,  texture  measures,  and  synthetic  structure  sizes)  was 
used  to  calculate  a  standard  deviation  of  approximately 
0.003.  In  the  CDF  difference  graphs  [labeled  (b)],  the  maxi¬ 
mum  difference  between  the  CDFs  is  equal  to  one.  Thus,  the 
lower  the  value,  the  better  the  agreement  between  the  clinical 
and  synthetic  textures.  The  CDF  difference  is  translated  into 
a  significance  level,  as  explained  previously,  on  the  vertical 
axis  at  the  right  of  the  graphs.  Values  of  the  CDF  difference 
corresponding  to  significance  levels  of  1%,  5%,  and  20%  are 
shown. 

It  can  be  seen  that  the  size  analysis  (represented  by 
circles)  shows  better  agreement  between  the  clinical  and  syn- 
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Fig.  7.  Summary  of  (a)  histogram  correlation  and  (b)  the  results  of  the  KS 
test,  for  subcutaneous  adipose  tissue.  The  data  are  presented  for  three  ranges 
of  size  of  synthetic  tissue  structures  and  for  the  three  texture  analysis  meth¬ 
ods. 


(b) 


Fig.  8.  Summary  of  (a)  histogram  correlation  and  (b)  the  results  of  the  KS 
test,  for  retromammary  adipose  tissue.  The  data  are  presented  for  three 
ranges  of  size  of  synthetic  tissue  structures  and  for  the  three  texture  analysis 
methods. 


thetic  texture,  than  the  other  two  features.  This  is  expected  as 
the  size  of  the  radiolucencies  is  related  to  the  size  of  the 
adipose  compartments,  while  the  other  features  are  more  sen¬ 
sitive  to  local,  small  scale  structure.  The  current  version  of 
our  model  does  not  include  fine,  local  tissue  detail. 

Results  for  the  retroareolar  glandular  tissue  regions  are 
shown  in  Fig.  9.  Most  of  the  feature  values  are  concentrated 
very  close  to  one  in  the  histogram  correlation  graph,  indicat¬ 
ing  good  agreement  between  the  simulation  and  the  real 
mammograms.  Also,  the  CDF  difference  for  the  retroareolar 
glandular  tissue  are  lower  that  for  the  other  tissue  regions. 
The  dense  glandular  region  shows  similarly  good  agreement 
(Fig.  10).  By  comparison,  both  glandular  tissue  regions  are 
simulated  better  than  the  adipose  regions  (Figs.  7  and  8). 

A  repeated  comparison  between  synthetic  and  real  mam¬ 
mograms  at  a  resolution  of  400  /tm/pixel  (not  plotted) 
showed  that  the  simulation  results  were  not  affected  signifi¬ 
cantly  by  the  change  of  resolution.  In  addition,  from  Figs.  7 
to  10  one  can  see  that  the  analyzed  features  are  sensitive  to 
the  size  of  simulated  anatomical  structures  and  that  the 
agreement  between  the  synthetic  and  clinical  mammograms 
depends  upon  our  selection  of  the  simulated  structure  size.  A 


partial  analysis  of  a  larger  set  of  synthetic  mammograms  was 
also  performed  and  no  significant  changes  in  the  comparison 
with  the  set  of  real  images  were  found. 

The  agreement  between  distributions  of  mean  texture  fea¬ 
tures  suggests  that  the  synthetic  images  sufficiently  well 
simulate  mean  features  of  the  clinical  texture.  Similarity  be¬ 
tween  the  averaged  histograms  of  real  and  synthetic  texture 
features  means  that  our  synthetic  images  can  simulate  the 
range  of  features  seen  over  a  large  group  of  mammograms, 
not  necessarily  matching  the  feature  distribution  of  any  par¬ 
ticular  mammogram. 

By  varying  the  parameters  which  control  the  sizes  of 
breast  tissue  model  elements,  we  were  able  to  match  the 
average  statistical  properties  of  clinical  mammograms  for  all 
tissue  types  except  the  retromammary  fat.  The  best  match  for 
clinical  mammogram  texture  was  achieved  for  the  simulated 
compartments  with  radii  of  4-13.3  mm  (“Medium"  and 
“Large,"  as  labeled  in  Table  I  in  the  accompanying  paper1) 
in  predominantly  adipose  tissue  region,  and  with  radii  of 
2.7-5.33  mm  (“Large")  and  1. 3-2.7  mm  (“Small")  in  the 
retroareolar  and  dense  FGT  region,  respectively.  These  pa- 
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Fig.  9.  Summary  of  (a)  histogram  correlation  and  (b)  the  results  of  the  KS 
test,  for  retroareolar  fibroglandular  tissue.  The  data  are  presented  for  three 
ranges  of  size  of  synthetic  tissue  structures  and  for  the  three  texture  analysis 
methods. 


(b) 


Fig.  10.  Summary  of  (a)  histogram  correlation  and  (b)  the  results  of  the 
KS  test,  for  dense  fibroglandular  tissue.  The  data  are  presented  for  three 
ranges  of  size  of  synthetic  tissue  structures  and  for  the  three  texture  analysis 
methods. 


rameters  were  chosen  because  of  high  histogram  correlation 
values  and  low  maximum  CDF  difference;  the  corresponding 
distribution  of  the  size  analysis  features  for  real  and  synthetic 
images  cannot  be  distinguished  at  the  5%  level.  A  difference 
between  the  retroareolar  and  dense  FGT  region  is  expected, 
since  the  retroareolar  region  contains  more  fat  clustered  in 
larger  compartments  than  in  dense  regions.  Further  under¬ 
standing  of  the  clinical  retromammary  adipose  tissue  struc¬ 
ture  is  needed  to  improve  the  simulation. 

IV.  CONCLUSIONS 

Evaluation  of  the  synthetic  mammograms  was  performed 
by  texture  analysis  and  comparison  with  normal  clinical 
mammograms  from  the  MIAS  database.  By  varying  the  dis¬ 
tribution  of  tissue  structures  in  the  model  we  have  been  able 
to  match  some  of  the  statistical  properties  of  clinical  mam¬ 
mograms.  Quantitatively,  the  synthetic  mammograms  have  a 
similar  distribution  of  the  values  averaged  over  a  large  num¬ 
ber  of  mammograms  for  several  texture  features,  namely  the 
average  size  of  image  objects,  the  texture  energy,  and  the 
fractal  dimension.  The  analysis  of  mammogram  object  size  is 


closely  related  to  the  size  analysis  of  the  medium  scale  phan¬ 
tom  elements,  simulating  adipose  compartments  in  the 
breast,  which  are  responsible  for  generation  of  the  synthetic 
mammogram  texture.  This  was  the  first  such  use  of  the  mor¬ 
phological  analysis  of  lucent  mammogram  regions,  repre¬ 
senting  the  projections  of  adipose  compartments.  Previous 
applications  of  morphology  focused  on  the  mammogram  re¬ 
gions  brighter  than  their  surroundings,  such  as  the  fibrous 
structures  and  microcalcifications.  The  texture  energy  and 
the  fractal  dimension  are  more  sensitive  to  the  local  variation 
of  pixel  intensities  due  to  the  small  scale  breast  tissue  struc¬ 
tures. 

Our  model  clearly  captures  the  coarse  tissue  structures  of 
the  breast  for  all  the  tissue  subregions  except  the  retromam¬ 
mary  fat.  This  exception  is  likely  due  to  the  fact  that  cur¬ 
rently  the  retromammary  adipose  tissue  is  simulated  in  the 
same  way  as  the  subcutaneous  tissue,  although  it  is  possible 
that  the  amount  of  fibroglandular  tissue  differs  in  these  two 
regions.  The  model  is  less  capable  of  capturing  the  small 
scale  structures  of  the  breast,  e.g.,  blood  and  lymph  vessels 
and  fine  tissue  detail,  which  affect  fine  texture  and  give  or- 
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ganized  structure  familiar  to  radiologists.  The  overly  geomet¬ 
ric  appearance  of  the  borders  between  the  AT  and  FGT  re¬ 
gions  in  the  synthetic  mammograms  can  be  improved  by 
small,  random  variations  in  the  position  of  the  borders  in  the 
compressed  tissue  model.  Differences  between  the  synthetic 
and  clinical  images  are  more  evident  for  texture  measures 
that  emphasize  smaller  spatial  scales,  in  agreement  with  the 
qualitative  visual  assessment.  We  expect  that  the  introduction 
of  detailed  tissue  structures  in  our  breast  model  will  enhance 
the  local  variations  of  synthetic  mammograms  and  the  varia¬ 
tions  in  feature  distribution  needed  to  better  match  clinical 
images. 
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APPENDIX:  BACKGROUND  ON  TEXTURE 
ANALYSIS 

1.  Mathematical  morphology 

Morphological  operations  are  performed  on  an  image, 
/(jc,y),  using  a  second  set  of  pixels,  S,  called  the  structuring 
element.  The  basic  morphological  operations,  are  defined 
by:13 

Dilation:  (/© S) (x,y )  =  ma x{f(x  + a;  \y + y ' )  |  (x  \y ' )  e  S}  > 

Erosion:  (fQS)(x,y)  =  Tmn{f(x-x' ,y-y')\(x\yr)  eS}, 

(Al) 

Opening:  /OS  =  (/0S)©S, 

Closing:  fS  =  (f®S)eS . 

2.  Texture  energy  analysis 

2D  convolution  masks  for  texture  energy  analysis  are  de¬ 
rived  using  ID  masks  with  different  number  of  zero- 
crossings,  designed  to  detect  different  texture  properties.  For 
example,  the  five-element  ID  masks  are:9 

Level:  L5  =  [l  4  6  4  1], 

Edge:  £5  =  [-l  -2  0  2  1], 

Spot:  S5  =  [- 1  0  2  0  -1], 

Wave:  W5=[- 1  2  0-2  1], 

Ripple:  R5=[i  -4  6  -4  1]. 

The  most  often  used  2D  masks  in  mammogram  analysis 
are  5X5  “lever'  and  “ripple"  masks,  obtained  by  the  prod¬ 
uct  of  the  corresponding  ID  masks  (R5R5=R5rR5  and 
L5L5=L5rL5): 
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3.  Fractal  analysis 

Image  fractal  dimension  is  usually  defined  using  the  area 
of  the  image  surface.  When  the  scale,  e,  is  increased  (which 
corresponds  to  decreasing  the  resolution)  the  area  of  a  fractal 
surface,  A  (e),  decreases.  The  fractal  dimension,  Z),  is  related 
to  the  slope  of  decreasing  area  on  a  log-log  plot,  as 

log  A(e)  =  log  const©  (2  -D)log  e.  (A3) 

There  are  several  algorithms  for  computing  fractal  dimension 
based  upon  box  counting,  image  power  spectrum,  or  iterated 
function  systems.26  Any  attempt  at  measuring  fractal  dimen¬ 
sion  must  deal  with  the  fact  that  self-similarity  of  real  images 
holds  only  over  a  limited  range  of  scales,  due  to  actual  struc¬ 
ture  and  limitations  of  the  imaging  process.27 
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Classification  of  Galactograms 
with  Ramification  Matrices: 

Preliminary  Results1 
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Rationale  and  Objectives.  The  poor  specificity  of  galactography,  the  imaging  modality  generally  indicated  in  cases  of 
^nnle  d  scWc  has  kd  o  a  large  number  of  biopsies  with  negative  results.  A  quantitative  scheme  for  classing  galac- 

<5^ 

ces)  and  the  correlation  of  the  values  of  the  matnx  elements  with  clinical  findings. 

Materials  and  Methods  The  ductal  trees  were  manually  segmented  for  25  galactographic  views  from  15  patients,  and 

matrices  of  the  pruned  trees.  ,  ..  ,  .  . 

n  matriv  oipment  Q2%  and  62%  of  NF  and  RF  cases  were  correctly  classified,  respectively 

|*;S"  *007) 'with  use  of  two  elements,  83%  and  77%  of  NF  and  RF  cases  were  correctly  classified,  but  this  result  was  not 
statistically  significant  (P  =  .108).  In  a  test  of  robustness,  an  analysis  of  pruned  trees  yielded  an  average  rool-mean-sq 
fractional  difference  of  9.7%  between  the  elements  of  the  original  and  the  R  matrix  average  over  a  prune 
Conclusion.  The  preliminaiy  analysis  indicates  that  it  may  be  possible  to  identify  cases  with  reported  galactographic  find¬ 
ings  by  using  R  matrices. 

Key  Words.  Breast,  ducts;  breast  neoplasms,  diagnosis;  galactography. 


Virtually  all  breast  cancers  derive  from  epithelial  tissue, 
with  90%  of  malignant  lesions  arising  in  the  ductal  epi¬ 
thelium  (1,  p  118).  Most  carcinomas  spread  initially  along 
the  lumen  of  the  ducts  or  lobules.  Evidence  of  this 
growth  pattern  is  seen  in  the  distribution  of  calcifications 
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associated  with  early  breast  cancer,  which  often  follows 
the  ductal  pathway  (2). 

Major  ducts  of  the  breast  extend  from  the  nipple  to¬ 
ward  the  chest  wall  in  a  branching  network  of  smaller 
and  smaller  ducts,  which  defines  a  draining  territory,  or 
lobe  (1,3-6).  The  adult  breast  contains  15-20  irregular 
lobes,  which  converge  to  the  nipple.  Each  lobe  is  drained 
by  its  own  major  duct.  Several  major  ducts  merge  to  form 
an  ampulla  (or  lactiferous  sinus),  a  dilated  segment  be¬ 
neath  the  nipple.  Branching  of  the  ducts  toward  the  chest 
wall  continues  until  a  duct  finally  ends  in  blunt  fingerhke 
ductules  formed  by  the  acini,  the  basic  glandular  secre¬ 
tory  units.  During  lactation  the  ductal  network  drains  milk 
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produced  in  the  acini.  Ductules  surrounded  by  specialized 
connective  tissue  are  called  lobules.  A  lobule  with  its  ter¬ 
minal  duct  is  known  as  the  terminal  ductal  lobular  unit. 

Some  cancers  are  revealed  by  nipple  discharge  alone, 
with  no  palpable  or  mammographic  lesions.  Indicated  in 
cases  of  nipple  discharge,  galactography  is  a  procedure 
for  imaging  the  contrast  material- enhanced  ductal  net¬ 
work  (7-9).  It  is  performed  by  carefully  identifying  the 
discharging  nipple  orifice,  introducing  a  blunt  needle,  and 
injecting  a  small  amount  of  radiographic  contrast  mate¬ 
rial.  Pre-  and  postcontrast  mammograms  are  obtained 
with  the  needle  in  place,  thereby  revealing  the  breast  lobe 
that  contains  the  discharging  duct.  Various  ductal  patterns 
(eg,  filling  defects,  ductal  ectasia)  can  be  recognized  from 
galactograms  (10).  Galactography  lacks  specificity  (1),  a 
situation  that  results  in  a  large  number  of  biopsies  with 
normal  or  benign  results.  A  quantitative  radiographic  clas¬ 
sification  scheme  for  galactograms  might  help  reduce  the 
percentage  of  biopsies  with  negative  results  and  the  re¬ 
lated  psychological  and  economic  effects. 

The  ductal  origin  of  breast  cancer  is  the  physiological 
basis  for  various  techniques  of  diagnosis  and  treatment. 
The  analysis  of  nipple  aspirate  fluid  (1 1)  has  been  investi¬ 
gated  for  early  cancer  detection.  Breast  cancer  risk  has 
been  estimated  by  analyzing  the  parenchymal  pattern  of 
projected  fibrous  and  ductal  structures  (12).  In  surgical 
assessment  of  nipple  discharge,  preoperative  staining  of 
the  ducts  can  minimize  the  amount  of  tissue  excised  (9). 
The  breast  ductal  network  has  been  modeled  by  tracing 
the  points  of  duct  entrance  to  and  exit  from  a  series  of 
subgross  histologic  slices  (13,14).  The  use  of  statistical 
analyses  of  ductal  networks  for  breast  modeling  also  has 
been  reported  (15,16). 

Describing  the  normal  ductal  network  is  difficult  be¬ 
cause  of  anatomic  variability  and  low  radiographic  con¬ 
trast.  Here,  we  report  a  method  of  analyzing  ductal  net¬ 
works  by  using  ramification  matrices  (R  matrices),  which 
describe  the  topologic  shape  of  a  treelike  structure  (17). 
The  R-matrix  elements  represent  branching  probabilities 
at  various  levels  of  a  tree  and  can  be  used  to  describe  a 
given  tree  or  to  generate  a  family  of  trees.  The  probabilis¬ 
tic  nature  of  R  matrices  makes  them  useful  for  generating 
many  individual  synthetic  trees  with  matching  statistical 
properties  of  branching.  For  this  reason,  R  matrices  have 
been  used  to  generate  ductal  networks  in  breast  modeling 
for  mammographic  simulation  (16). 

The  morphology  of  the  ductal  network  reflects  the 
state  of  breast  development  and  the  healthy  or  pathologic 
state  of  the  breast  tissue,  as  has  been  shown  by  analyses 


of  the  branching  of  murine  mammary  ducts  in  different 
phases  of  gland  development  (18,19)  and  by  studies  of 
epithelial  cellular  organization  under  the  influence  of  hor¬ 
mones,  growth  factors,  and  carcinogens  (20-22).  We  hy¬ 
pothesize  that  diseases  of  the  breast,  demonstrated 
through  alterations  of  the  normal  ductal  anatomy,  can  be 
quantified  and  classified  from  galactograms.  In  the  present 
study,  we  manually  segmented  ductal  trees  from  25  galac¬ 
tograms,  calculated  the  R  matrices,  and  computed  the 
probability  of  correct  classification  for  a  combination  of 
R-matrix  elements.  We  also  tested  the  robustness  of  the 
R-matrix  representation  by  analyzing  pruned  trees. 


MATERIALS  AND  METHODS 


Clinical  galactograms  for  this  analysis  were  obtained 
retrospectively  from  15  patients  with  a  mean  age  of  49.2 
years  (range,  29-75  years).  These  patients  were  selected 
from  a  group  of  41  who  had  undergone  galactography  at 
the  Thomas  Jefferson  University  Breast  Imaging  Center, 
Philadelphia,  Pa,  during  the  6'/2-year  period  from  June 
1994  through  January  2001.  Galactograms  from  17  of  the 
41  patients  were  unavailable  because  they  had  been  re¬ 
turned  to  the  patients  or  the  primary  health  care  institu¬ 
tions,  and  galactograms  from  another  nine  patients  were 
not  used,  because  of  obstruction  or  poor  image  quality 
(ie,  the  complete  ductal  tree  could  not  be  segmented).  Of 
the  15  patients  whose  cases  were  analyzed,  eight  (mean 
age,  44.2  years;  range,  29-74  years)  had  no  reported  ga- 
lactographic  findings  (hereafter  NF),  and  seven  (mean 
age,  54.8  years;  range,  43-75  years)  had  galactographic 
findings  of  ductal  ectasia,  cysts,  or  papilloma  (hereafter 
RF).  Twenty-five  galactographic  views  of  the  15  patients 
were  analyzed  (16  craniocaudal  [CC]  and  nine  mediolat- 
eral  or  mediolateral  oblique  [hereafter  denoted  in  combi¬ 
nation  as  ML/MLO]),  of  which  12  views  (eight  CC  and 
four  ML/MLO)  were  from  NF  cases  and  13  (eight  CC 
and  five  ML/MLO)  were  from  RF  cases.  There  were  no 
reported  findings  of  malignancy  from  these  25  galacto¬ 
grams.  Furthermore,  mammographic  follow-up  data  were 
available  for  eight  of  the  15  patients  for  an  average  pe¬ 
riod  of  4.75  years  (range,  3.5-6.0  years),  and  no  malig¬ 
nancies  were  reported.  Patients’  ages,  available  galacto¬ 
graphic  views,  symptoms,  and  reported  galactographic 
findings  are  listed  in  the  Table. 

To  reconstruct  the  ductal  topology,  each  branch  in  the 
ductal  network  was  drawn  by  hand  on  a  sheet  of  tracing 
paper  placed  over  the  galactogram  displayed  on  a  light 
box  (Figure  1).  The  points  where  ducts  branched  were 
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Age  (y) 

View 

Symptom 

Galactographic  Finding 

29 

LCC 

Greenish  discharge 

None 

30 

LML,  LCC 

Clear/yellow  discharge 

None 

32 

LML(Mag),  LCC 

Bloody  discharge 

None 

36 

LCC 

Bloody/yellow  discharge 

None 

43 

RML,  RCC 

Greenish  discharge 

Cysts 

44 

RCC 

Bloody  discharge 

Cysts,  ductal  ectasia 

45 

RCC,  LML 

Bloody/milky  discharge 

None 

45 

LCC 

Dilated  ducts  on  mammogram  and  US  image 

None 

47 

RMLO,  RCC 

Greenish  discharge,  mastitis 

Cysts,  ductal  ectasia 

RML,  RCC 

50 

RML,  RCC 

Greenish  discharge 

Cysts 

55 

RCC,  RMLO 

Clear  discharge 

Cysts 

63 

LMLO,  LCC 

Darkish  discharge 

None 

70 

RCC 

Bloody  discharge 

Ductal  ectasia 

74 

LCC 

Clear/bloody  discharge 

None 

75 

RCC 

Bloody  discharge 

Papilloma 

Note.— LCC  =  left  craniocaudal,  LML  =  left  mediolateral,  RML  =  right  mediolateral,  RCC  =  right  craniocaudal,  RMLO  -  right  medio- 
lateral  oblique,  LMLO  =  left  mediolateral  oblique,  and  Mag  =  magnified. 


distinguished  from  the  points  of  overlap  by  the  fact  that 
the  latter  are  galactographically  brighter,  due  to  superpo¬ 
sition  of  the  x-ray  attenuation.  Large  ducts  were  recon¬ 
structed  by  connecting  the  marked  points.  In  each  of  the 
segmented  ductal  trees,  the  root,  internal  and  terminal 
nodes,  and  branches  were  labeled  and  the  R-matrix  ele¬ 
ments  were  computed  as  described  by  Viennot  et  al  (17). 
Simplified,  the  algorithm  is  given  as  follows:  ( a )  all  ter¬ 
minal  branches  have  label  1,  (b)  a  “parent”  branch  whose 
“children”  have  labels  i  and  j  will  be  labeled  by  max(i,  j) 
if  /  +  j  or  by  (i  +  1)  if  i  =  j,  and  (c)  the  labeling  proce¬ 
dure  continues  until  the  root  branch  is  reached  whose  la¬ 
bel  s  is  called  the  Strahler  number  of  the  tree  structure. 

The  R  matrix  of  a  tree  with  Strahler  number  s  is  a 
lower  triangular  matrix,  defined  as 

=  [/•*,,  =  bkJlak,  k  G  (2,  s),j<=(  1 ,  *)],  (1) 

where  ak  is  equal  to  the  number  of  branches  with  label  k 
(17).  For  j  <  k,  bkJ  is  the  number  of  pairs  of  branches 
with  labels  k  and  j,  while  for  j  =  k,  bkJ  is  the  number  of 
pairs  of  branches  both  labeled  k  -  1,  descending  from  a 
node.  Therefore,  rkj  =  bkJak  =  p(bkj\ak)  is  the  probability 
that  a  branch  with  label  k  will  bifurcate  into  branches 
with  the  appropriate  labels.  Figure  lc  shows  the  numeri¬ 
cally  labeled  branches  of  the  segmented  ductal  tree  from 
Figure  la.  In  Figure  lc,  there  are  14  branches  with  label 
2,  six  of  which  bifurcate  into  pairs  of  branches  with  la¬ 


bels  1  and  2,  corresponding  to  the  probability  of  r2>i  = 
6/14  =  0.43.  The  other  eight  branches  with  label  2  bifur¬ 
cate  into  pairs  of  branches  both  with  label  1 ,  correspond¬ 
ing  to  the  probability  of  r2, 2  =  8/14  =  0.57.  In  a  similar 
manner,  all  the  elements  of  the  R  matrix  shown  in  Figure 
Id  were  computed  from  the  galactogram  shown  in  Figure 
la.  We  computed  R  matrices  with  nine  elements  and 
Strahler  number  s  =  4,  corresponding  to  a  root  branch 
with  label  4. 

A  linear  Bayesian  decision  rule  was  used  to  classify 
the  galactographic  findings.  In  the  design  of  the  Bayesian 
classifier,  we  assumed  that  the  data  were  normally  distrib¬ 
uted  and  that  the  population  standard  deviation  was  the 
same  for  both  classes.  These  assumptions  reduce  the  risk 
of  overspecification  due  to  the  small  sample  size.  We  es¬ 
timated  the  standard  error  of  the  classification  results  by 
means  of  a  leave-one-out  (jackknife)  method,  in  which 
the  percentages  correctly  classified  are  calculated  for  sub¬ 
sets  of  the  data  formed  by  leaving  each  sample  out  in 
turn.  The  standard  error  of  the  classification  results  is 
then  the  standard  deviation  of  the  values  calculated  for 
the  jackknife  subsamples  (23). 

A  preliminary  test  of  the  robustness  of  the  R-matrix 
approach  was  performed  by  analyzing  pruned  ductal  trees 
derived  from  an  original  tree  after  removal  of  a  single 
terminal  branch.  This  test  was  proposed  because  of  the 
observed  ambiguity  in  tracing  terminal  branches,  which 
can  be  easily  overlooked  due  to  their  small  size  and  low 
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a. 

b. 

C. 

_  -i 

1 

7*2,1  7*2,2 

0.43  0.57  . 

R  = 

7*3.1  7*3,2  7*3,3  * 

= 

0  0.33  0.67  . 

r.i,i  7*4,2  7*4,3  7*4,4 

0  0.75  0  0.25 

d. 


Figure  1 .  Segmentation  of  a  ductal  tree, 
showing  (a)  part  of  a  galactogram  with  a  con¬ 
trast-enhanced  ductal  network,  (b)  the  manu¬ 
ally  traced  network  of  larger  ducts  from  the 
contrast- enhanced  portion  of  the  galacto¬ 
gram,  (c)  numeric  labeling  of  branches  in  the 
ductal  network,  and  (d)  the  R  matrix  com¬ 
puted  from  the  branching  pattern.  The  dots, 
triangles,  and  squares  denote  branching 
points  of  different  levels  of  the  tree. 


contrast.  Consequently,  some  of  the  terminal  branches 
might  have  been  omitted,  resulting  in  a  pruned  version  of 
the  original  tree.  The  R  matrix  of  a  pruned  tree  was  com¬ 
puted  and  the  procedure  was  repeated  for  each  of  the  ter¬ 
minal  branches  in  the  original  tree.  A  comparison  was 
performed  by  computing  the  root-mean-square  fractional 
difference,  as  follows: 


D  = 


r  _  -prun" 
r  k.j  ~  k,j 


n,j 


(2) 


where  rkJ  and  are  the  elements  of  the  original  R  ma¬ 
trix  and  the  matrix  averaged  over  all  pruned  versions  of 
the  original  tree,  respectively.  Nei  is  the  number  of  non¬ 
zero  R-matrix  elements  (A^i  —  9  for  s  =  4,  from  Eq  [1]). 


RESULTS 


The  values  of  the  R-matrix  elements  calculated  from 
clinical  galactograms  and  the  averages  over  all  NF  and  all 


RF  cases  are  plotted' in  Figure  2.  Some  matrix  elements 
show  a  noticeable  difference  in  mean  value  between  the 
two  classes  (eg,  r3>2  and  r3t3),  suggesting  the  possibility  of 
classifying  galactographic  findings  on  the  basis  of  R  ma¬ 
trices.  We  evaluated  a  classification  scheme  that  used  ei¬ 
ther  a  single  matrix  element  or  a  pair  of  them  and  a  linear 
Bayesian  decision  rule  in  a  leave-one-out  fashion.  The 
classification  performance  of  a  single  R-matrix  element 
was  best  (in  terms  of  the  sum  of  the  correctly  identified 
fractions  in  the  two  groups)  for  r3>3,  which  correctly  clas¬ 
sified  92%  ±  2  of  NF  and  62%  ±  3  of  RF  cases.  In  a 
combination  of  two  elements,  r3>2  and  r3>3  correctly  classi¬ 
fied  83%  ±  4  of  NF  cases  and  77%  ±  4  of  RF  cases. 
Figure  3  illustrates  the  classification  results. 

To  estimate  the  statistical  significance  of  the  classifica¬ 
tion  results,  we  performed  an  additional  experiment  by 
using  random  numbers  as  descriptive  parameters.  Six 
numbers,  corresponding  to  the  linearly  independent  ele¬ 
ments  of  the  R  matrix  with  nine  elements,  were  randomly 
generated  for  each  of  the  25  galactograms.  This  procedure 
was  repeated  for  10,000  trials,  and  in  each  trial  the  ran- 
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all  NF  cases,  and  light  symbols  represent  the  individual  cases. 
Error  bars  correspond  to  the  sample  standard  deviations. 


domly  generated  data  were  used  to  classify  galactograms 
by  means  of  either  one  or  two  matrix  elements,  as  was 
done  for  the  real  data.  In  0.7%  of  trials,  one  of  the  pa¬ 
rameters  showed  a  classification  power  of  at  least  92% 
for  one  class  and  at  least  62%  for  the  other— comparable 
with  the  performance  of  element  r3i 3  on  the  clinical  galac¬ 
tograms.  Thus,  if  the  classification  result  observed  for  the 
clinical  galactograms  is  by  chance  alone,  it  would  be  en¬ 
countered  in  just  0.7%  of  repeated  experiments,  and  the 
statistical  significance  for  classification  with  r 3,3  is  P  = 
.007.  The  classification  results  achieved  with  the  pair  of 
elements  r3i2  and  r3>3  was  not  statistically  significant  (P  = 
.108). 

The  analysis  of  the  pruned  versions  of  the  clinical  ga¬ 
lactograms  showed  an  average  root-mean-square  frac¬ 
tional  difference  between  the  original  R-matrix  elements 
and  the  elements  averaged  over  all  pruned  trees  of  9.7% 
for  all  NF  and  RF  galactograms  (6.5%  for  the  NF  cases 
and  12.6%  for  the  RF  cases).  A  similar  analysis,  which 
considered  1,000  synthetically  generated  random  binary 


Figure  3.  Classification  of  galactographic  findings  based  on  ele¬ 
ments  of  the  R  matrix  and  a  Bayesian  linear  decision  rule.  Bold 
symbols  (O  =  RF,  *  =  NF)  represent  element  values  averaged 
over  all  RF  and  NF  cases,  and  light  symbols  represent  the  individ¬ 
ual  cases.  Short  dashes  indicate  Bayesian  decision  line  for  r3, 3, 
and  long  dashes  indicate  Bayesian  decision  line  for  r3,2,  r3i3. 

trees  pruned  in  the  same  fashion,  yielded  an  average  root- 
mean-square  fractional  difference  of  6.8%  (24). 


DISCUSSION 


We  investigated  the  branching  structure  of  ducts  visu¬ 
alized  on  galactograms  and  evaluated  the  use  of  R  matri¬ 
ces  to  classify  galactographic  findings.  As  an  illustration 
of  the  classification  results,  Figure  4  shows  two  galacto¬ 
grams  from  the  sets  of  NF  and  RF  cases  used  in  the 
study.  The  NF  galactogram  in  Figure  4a  corresponds  to 
the  values  r3,2  =  0.5  and  r„  =  0.19.  The  RF  galactogram 
in  Figure  4b  corresponds  to  the  values  r3  2  =  0.33  and 
r3,3  =  0.67. 

Results  obtained  in  this  study  show  that  element  r3,2  is 
approximately  50%  smaller  and  element  r 3  3  is  approxi¬ 
mately  50%  larger  when  averaged  over  all  RF  images, 
relative  to  the  average  over  NF  images:  <r3>2>  =  0.20 
whereas  <f3I,2>  =  0.38,  and  <r'3  3>  =  0.52  whereas 
<r3,3>  =  0.33.  These  differences  in  the  matrix  elements 
may  be  due  to  long-term  processes  that  alter  the  ductal 
branching  pattern.  Another  explanation  may  be  cutoff  or 
dilation  of  the  ducts,  which  affects  the  filling  of  the  duc¬ 
tal  network  with  contrast  agent  and  the  visibility  of  the 
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Figure  4.  Two  examples  of  galactograms 
that  have  been  correctly  classified  by  means 
of  R  matrices,  (a)  Galactogram  with  no  re¬ 
ported  findings  (patient  age,  45  years;  right 
CC  view;  r3i2  =  0.5  and  r3i3  =  0.19).  (Large 
bright  regions  seen  in  this  galactogram  are 
due  to  extravasation,  which  did  not  affect  the 
segmentation  of  the  ductal  tree.)  (b)  Galacto¬ 
gram  with  a  reported  finding  of  cysts  (patient 
age,  55  years;  right  CC  view;  r3i2  =  0.33  and 
r 3f3  —  0.67). 


a. 


b. 


smaller  ducts,  thereby  altering  the  calculated  R  matrix. 
Further  research  is  needed  to  investigate  the  relationship 
between  galactographic  changes  and  R-matrix  element 
values. 

The  root-mean-square  fractional  difference,  estimated 
from  analysis  of  the  pruned  ductal  trees,  is  on  average 
9.7%,  significantly  smaller  than  the  difference  between 
the  average  values  of  matrix  elements  r3>2  and  r3>3.  More¬ 
over,  the  sample  standard  deviation  of  these  elements  is 
significantly  larger  than  the  root-mean-square  difference 
measured  from  pruning.  The  means  and  standard  devia¬ 
tions  of  these  elements  are  shown  in  Figure  2.  These 
findings  substantiate  the  robustness  of  the  R-matrix  repre¬ 
sentation  of  ductal  networks,  as  the  observed  variation  in 
the  galactograms  cannot  be  explained  by  the  accuracy  of 
the  segmentation. 

Evaluation  with  the  simulated  data  showed  that  clas¬ 
sification  based  on  a  single  element,  r3>3,  performs  at  a 
statistically  significant  level  (P  =  .007).  The  perfor¬ 
mance  of  classification  based  on  a  pair  of  elements,  r32 
and  r3  3,  was  not  significant  (P  =  .108).  This  reduction 
in  significance  is  understandable  given  that  the  use  of 
two  matrix  elements  increases  the  probability  that  some 


pair  will  give  a  spuriously  efficacious  classification 
and  given  that  the  classifier  design  assigned  equal  im¬ 
portance  to  each  matrix  element  to  avoid  overspecifica¬ 
tion. 

Several  confounding  factors  should  be  considered  in 
the  interpretation  of  these  results.  Patient  age  distributions 
differ  for  the  NF  and  RF  cases  (mean  age,  44.2  years  for 
NF  vs  54.8  years  for  RF).  We  tested  and  found  essen¬ 
tially  no  correlation  between  the  age  and  the  matrix  ele¬ 
ment  values,  with  Pearson  correlation  coefficient  values 
of  0.18  and  0.01  for  the  elements  r3>2  and  r3  3,  respec¬ 
tively.  Another  possible  influence  could  be  our  combina¬ 
tion  of  galactograms  obtained  with  all  views,  CC  and 
ML/MLO.  In  an  ideal  case,  with  all  the  branches  visible 
and  with  perfect  segmentation,  the  reconstructed  ductal 
trees  would  be  the  same  on  both  CC  and  ML/MLO 
views,  as  would  the  corresponding  R  matrices,  which  was 
not  the  case  in  this  experiment.  We  compared  the  maxi¬ 
mum  difference  between  the  mean  element  values  com¬ 
puted  over  normal  and  benign  galactograms  with  the 
maximum  difference  between  the  mean  element  values 
computed  over  CC  and  ML/MLO  views  and  found  the 
former  difference  approximately  two  times  greater  than 
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the  latter.  The  most  limiting  factor  in  our  study,  however, 
was  the  small  sample. 

In  summary,  our  preliminary  analysis  indicates  that  R 
matrices  may  be  used  to  identify  cases  with  reported  ga- 
lactographic  findings.  The  performance  of  classification 
with  a  single  matrix  element  was  statistically  significant 
in  a  set  of  eight  NF  cases  (12  views)  and  seven  RF  cases 
(13  views)  with  findings  of  ductal  ectasia,  cysts,  or  papil¬ 
loma.  Further  experiments  should  use  more  galactograms 
and  a  more  sophisticated  decision  rule  and  should  include 
malignant  cases. 
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I.  INTRODUCTION 


A  significant  limitation  of  projection  radiography  is  that  relevant  findings  are  often  ob¬ 
scured  or  mimicked  by  the  x-ray  shadows  of  other  anatomical  structures  (summation  ar¬ 
tifacts).  By  comparison,  stereoradiography  allows  superimposed  structures  to  be  spatially 
separated,  reducing  the  confounding  effects  of  overlap.  Stereoradiography  was  first  proposed 
by  Thomson  in  18961  and  continued  to  be  used  for  certain  radiographic  procedures  until  the 
1980’s.  The  first  recorded  use  of  stereomammography  occurred  in  1930. 2  The  advantages 
and  disadvantages  of  stereoradiographic  techniques  are  discussed  in  several  radiographic  text 
books  and  articles,  e.g.  Curry  et  al..3  Improvements  in  image  quality  of  conventional  film- 
screen  radiographic  systems  combined  with  the  added  effort^  cost  and  patient  dose  led  to  a 
decline  in  the  use  of  stereoradiography  and  stereomammography.3  Development  of  digital  x- 
ray  detectors  and  softcopy  reading  has  resulted  in  a  renewal  of  interest  in  stereoradiography 
as  a  potentially  useful  method  of  reducing  summation  artifacts. 

Dose  in  radiography  is  constrained  by  the  imaging  task,  the  x-ray  quantum  mottle,  and 
the  detector  specifications.  Quantum  mottle  is  caused  by  the  fluctuation  in  the  number 
of  x-ray  quanta  produced  and  transmitted  through  the  tissue.  For  the  same  detector,  the 
fewer  the  x-ray  quanta,  the  greater  the  relative  fluctuation,  and  the  lower  the  signal-to-noise 
ratio  of  the  object  to  be  detected.  An  example  of  detector  limitations  is  the  characteristic 
(“H  and  D”)  curve,  which  determines  a  minimal  useful  exposure.  The  issues  of  dose  in 
stereoradiography  are  either  not  mentioned  in  textbooks4-6  or  it  is  suggested  that  a  stere¬ 
ographic  technique  requires  a  total  of  twice  the  exposure  of  a  single-projection  technique.7 
Signal  detection  theory  suggests  that,  under  quantum  limited  conditions,  an  ideal  observer 
would  be  able  to  combine  information  from  the  two  stereoscopic  views  and  thereby  suppress 
quantum  mottle.  Thus,  under  this  model,  two  stereoscopic  views,  each  acquired  with  one 
half  of  the  dose  of  a  single-projection  radiograph,  would  provide  the  same  level  of  resilience 
to  quantum  noise  as  a  single  projection  acquired  with  the  same  total  dose,  but  provide  the 
inherent  benefits  of  stereoscopic  viewing. 
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Human  observers  perform  certain  visual  tasks  more  efficaciously  when  fusing  information 
from  both  eyes.8,9  In  this  paper  we  are  concerned  with  the  tasks  of  contrast  perception  and 
noise  suppression,  which  are  complementary  to  the  well  described  advantages  of  stereoscopy 
in  depth  perception.  For  example,  threshold  measurements  of  sinusoidal  patterns  have  shown 
a  decrease  in  threshold  (i.e.,  an  increase  in  sensitivity)  of  a  factor  of  \[2  when  viewing  is 
performed  binocularly  as  opposed  to  monocularly.10  This  increase  in  sensitivity  is  precisely 
as  expected  for  an  ideal  observer  combining  the  signal  from  each  eye.  The  quantum  mottle  in 
x-ray  projections  is  different  from  the  sources  of  noise  in  these  experiments;  in  a  radiograph, 
a  given  realization  of  the  noise  is  permanently  recorded  at  the  time  of  acquisition  and  it  is 
presented  to  each  eye  as  a  fixed  pattern.  For  example,  we  note  that  flat  random-noise  fields 
presented  stereoscopically  cause  lustre11  (i.e.,  the  image  to  appear  to  shimmer)  due  to  an 
attempt  by  the  visual  system  to  stereoscopically  fuse  the  random  bright  and  dark  variations. 

The  question  of  threshold  detection  of  objects  stereoradiographically  is  different  from 
detection  in  other  stereoscopic  settings.  To  address  this  difference  and  examine  the  role  of 
dose  in  stereoscopy  we  performed  an  observer  study.  Observer  studies  of  stereoradiography 
have  been  reported  in  the  literature  previously.  Kundel  et  al.9  compared  stereoscopic  acu¬ 
ity  and  effects  of  monocular  depth  cues  under  direct  vision  and  in  stereofluoroscopy,  while 
Goodsitt,  Chan,  and  Hadjiiski12  studied  depth  perception  in  stereomammography.  Berkson 
et  al.6  compared  the  number  of  false  negatives/positives  in  mono-  and  stereoscopic  chest 
radiographs  and  Hsu  et  al.13  studied  detection  of  simulated  abnormalities  in  stereomam¬ 
mography. 

Our  observer  study  was  based  upon  the  hypothesis  that  in  the  merging  of  left  and  right 
eye  images  by  the  human  visual  system,  the  dose  needed  for  each  of  the  two  stereoradio¬ 
graphic  images  is  equal  to  one  half  of  the  dose  for  a  single  x-ray  image  viewed  monoscopically. 
Thus,  by  corollary,  we  hypothesize  that  viewing  a  stereoscopic  image  will  result  in  an  ef¬ 
fective  increase  in  SNR  by  a  factor  of  y/2.  A  series  of  contrast-detail  (C-d)  experiments 
were  performed  with  stereoradiographic  images  acquired  over  a  range  of  exposures.  The 
C-d  experiment  belongs  to  a  class  of  threshold- visibility  psychophysical  measurements.14  In 


3 


the  study,  observers  attempted  to  detect  details  in  a  C-d  phantom.  The  observers  scored 
the  phantom  monoscopically  (both  eyes  seeing  the  same  radiographic  projection)  and  stereo- 
scopically  (each  eye  seeing  a  separately  acquired  projection).  The  focus  of  the  study  was 
on  the  issue  of  whether  the  human  observer  could  combine  the  stereoscopic  images  in  a 
manner  so  as  to  reduce  the  apparent  quantum  mottle.  The  two  projections  in  each  stere¬ 
oradiographic  pair  were,  therefore,  acquired  separately,  although  with  the  same  geometry. 
No  depth  information  was  encoded,  i.e.,  all  objects  were  positioned  in  the  imaging  plane, 
to  avoid  the  confounding  issues  of  the  variable  ability  of  observers  to  fuse  stereo  pairs.  The 
observers’  performances  were  evaluated  in  terms  of  the  difference  in  number  of  details  seen 
stereo-  vs.  monoscopically,  effective  stereoscopic  dose  (equal  to  the  dose  needed  for  a  mono- 
scopic  viewing  with  the  same  performance  as  for  the  stereoscopy) ,  and  the  effective  decrease 
in  the  threshold  SNR  values. 

II.  MATERIALS  AND  METHODS 
A.  Image  acquisition 

Radiographic  images  of  an  RMI-180  mammographic  C-d  phantom  (Gammex  RMI,  Mid¬ 
dleton,  WI)  were  used  in  the  observer  study.  The  phantom  has  90  objects  aligned  in  nine 
columns.  Adjacent  objects  within  the  same  column  differ  in  diameter  by  a  factor  of  \/2, 
with  a  range  of  7.07-0.312  mm.  The  thickness  of  the  objects  in  adjacent  columns  also  differ 
by  a  factor  of  %/2,  with  a  range  of  1.0-0.062  mm.  The  phantom  projections  were  acquired 
with  a  DirectRay  flat  panel  digital  x-ray  detector  (Hologic/Direct  Radiography,  Newark, 
DE),15  without  a  grid.  Observers  were  presented  with  a  subregion  of  650  x  810  pixels  with 
the  original  resolution  of  147  fim  pixels,  including  all  the  phantom  details. 

Imaging  was  performed  at  60  kVp  and  at  six  mAs  stations  in  the  range  of  2-100  mAs, 
using  a  general  radiography  Bi-150  30/50  x-ray  tube  (Siemens,  Munich,  Germany)  with 
tungsten  target,  a  measured  HVL  of  1.34  mm  A1  (at  60  kVp),  and  a  Heliophos  5S  generator 
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(Siemens,  Munich,  Germany).  Additional  filtration  of  6cm  of  lucite  was  used  to  simulate 
tissue  and  appropriately  harden  the  beam.  The  block  of  lucite  was  positioned  near  the  x-ray 
focus  in  order  to  reduce  scatter  and  provide  for  uniform  beam  filtration.  The  phantom  was 
imaged  in  contact  with  the  detector  (112  cm  from  the  focal  spot)  and  was  placed  within 
a  wide  lucite  frame  of  equal  thickness,  providing  uniform  scatter  throughout  the  whole 
phantom  area.  At  each  mAs  station,  we  acquired  five  images  of  the  phantom,  as  well  as  ten 
dark  field  (x-rays  off)  and  ten  bright  field  (x-rays  on,  with  the  phantom  removed)  images. 
The  bright  and  dark  fields  were  used  for  correcting  the  gain  and  offset  variations  of  the 
individual  detector  elements  (dels). 


B.  Image  processing 

The  raw  acquired  phantom  images  were  corrected  to  compensate  for  the  effects  of  detector 
nonuniformity  and  adjusted  for  the  non-linearity  of  human  contrast  sensitivity. 


1.  Detector  nonuniformity 


A  standard  procedure  for  reducing  the  pixel  variations  due  to  the  gain  and  offset  vari¬ 
ations  of  dels  is  to  apply  a  dark  and  bright  field  correction.  Corrected  pixel  values  were 
computed  as: 


E{Pjj} 

E{Bitj}  -  E{Dhi) 


+  I 


new  ) 


(i) 


where  /  and  Icorr  are  the  original  and  corrected  pixel  values,  respectively;  E{D}  and  E{B} 
are  the  averages  of  the  dark  and  the  bright  fields,  respectively;  and,  K  and  7new  are  param¬ 
eters  transforming  the  range  (contrast)  and  the  mean  image  pixel  values.  Subscripts  i  and 
j  denote  the  position  of  the  pixel  in  the  image  array. 

The  existence  of  malfunctioning  dels,  whose  corresponding  pixel  values  are  not  propor¬ 
tional  to  the  incident  x-ray  flux  and  differ  significantly  from  their  neighboring  pixels,  was 
observed  and  the  corresponding  pixel  values  corrected.  Recently,  several  researchers  have 
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investigated  this  problem.  Aach  and  Metzler16  proposed  image  deconvolution  using  the  spec¬ 
tral  analysis  of  the  defect  image  and  Tang  et  a/.17  used  wavelet  analysis  for  identification 
and  interpolation  of  flat  panel  images  used  for  cone  beam  CT. 

In  our  case,  the  average  value  and  the  variance  of  all  the  dark  fields  were  computed  for 
each  del,  and  the  same  procedure  was  repeated  for  all  of  the  bright  fields.  First,  dels  for 
which  the  digital  value  was  constant  accross  all  of  the  dark-  and  bright-  field  images  were 
identified  as  non-functioning.  The  pixels  corresponding  to  the  non-functioning  dels  were 
replaced  by  the  average  of  their  immediate  neighbors.  Second,  we  computed  the  spatial 
variance  within  3x3  pixel  neighborhoods  of  the  phantom  images,  and  averaged  the  spatial 
variance  values  over  a  small  and  large  neighborhood  (of  size  3  and  5  pixels  in  diameter, 
respectively).  The  averaged  variance  values  were  then  tested  against  the  following  criterion: 


<  O’  >4771 

<  O  >lg 


>  c, 


(2) 


where  <  a  >sm  and  <  o  >ig  are  the  spatial  averages  of  the  values  of  the  pixel  standard 
deviations,  computed  over  the  small  or  large  neighborhood,  respectively,  and  C  is  the  crite¬ 
rion  value.  Pixels  satisfying  the  criterion  were  replaced  by  the  average  of  their  immediate 
neighbors.  Many  of  the  affected  dels  were  grouped  along  several  vertical  and  horizontal 
lines.  The  value  of  parameter  C  that  would  correct  all  of  the  pixels  along  these  lines  would 
also  produce  an  undesirably  high  number  of  replaced  pixels  randomly  distributed  over  the 
image:  Replacing  too  many  correct  pixels  might  change  the  statistics  of  the  noise  and  affect 
the  results  of  the  study.  Therefore,  we  identified  the  lines  of  malfunctioning  dels  in  images 
and  directly  applied  the  correction.  For  the  rest  of  the  image,  the  value  of  C  =  2.1,  which 
produced  only  a  small  number  of  undesired  replacements,  was  selected. 


2.  Non-linear  human  contrast  sensitivity 

Human  contrast  sensitivity  is  a  non-linear  function  of  luminance.  Relative  changes  in 
luminance  are  more  easily  noticed  in  bright  image  areas  than  in  dark  areas.18  Perceptual 
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linearization  is  an  image  transformation  that  adjusts  the  brightness  so  that  equal  changes 
in  the  quantity  being  displayed  (here,  x-ray  fluence  at  the  detector  surface)  •will  be  equally 
perceived.  The  contrast  sensitivity  of  the  human  visual  system  is  approximated  by  Barten’s 

I  ' 

model.19  The  model  was  derived  for  a  standard  target  of  a  2  deg  x  2  deg  square  filled  with  a 
horizontal  or  vertical  sinusoidal  modulation  of  4  cycles/deg,  placed  in  a  uniform  background 
of  the  mean  target  luminance.18  The  threshold  modulation  at  which  the  target  becomes  just 
visible  to  the  average  observer  defines  the  just-noticeable  difference,  JND,  at  the  luminance 
value  of  the  background.  The  interpolation  of  the  luminance  levels  corresponding  to  1023 
JNDs  is  given  by  the  grayscale  standard  display  function:18 


log10  L(Na) 


Zlobin  NAf 
l  +  ZUk(lnNAy 


(3) 


where  NA  and  L  are  the  JND  index  (i.e.,  the  number  of  just-noticeable  differences)  and  the 
corresponding  luminance  value,  respectively.  Coefficients  a;  and  bi  of  the  rational  polynomial 
interpolation  are  given  in  the  DICOM  standard.18  We  have  incorporated  the  grayscale  stan¬ 
dard  display  function  into  the  transformation  from  the  flat-panel  output  12-bit  pixel  value 
to  the  monitor  input  8-bit  digital-driver  value.  This  transformation  consists  of  three  parts.20 
This  decomposition  offers  flexibility  to  perform  observer  studies  using  different  monitors  or 
at  different  overall  brightness  or  contrast  levels. 

In  the  first  step,  the  monitor  characteristic,  relating  luminance  levels  to  the  8-bit  digital 
driver  levels,  was  interpolated  from  photometric  measurements  made  by  a  TEK  Lumacolor 
J17  photometer  (Tektronix,  Beaverton,  OR).  The  luminance  values  were  measured  in  the 
center  and  at  the  periphery  of  the  monitor.  The  average  values  at  24  digital  driver  levels 
were  used  for  interpolation.  A  fourth  order  polynomial  fit  was  used  for  interpolation  with 
maximum  error  of  1.68  digital  driver  levels  (0.6%).  Second,  the  transformation  from  the 
luminance  values  to  the  number  of  JNDs  was  approximated  by  the  grayscale  standard  display 
function,  given  by  Eq.  3.  Third,  the  linear  transformation  from  12-bit  flat  panel  output  to 
the  JND  index  was  used  to  adjust  the  common  overall  brightness  and  contrast  of  all  the 
phantom  images  included  in  the  study.  The  overall  brightness  was  selected  in  a  small 
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preliminary  observer  study  as  the  one  which  gave  the  largest  increase  in  number  of  details 
seen  stereo-  vs.  monoscopically.  Based  on  that  study,  the  mid-level  between  the  largest, 
highest  contrast  element  and  its  annular  neighborhood  was  transformed  to  a  JND  index  of 
350,  corresponding  to  a  digital  driver  level  of  90. 

We  have  used  two  strategies  for  displaying  the  phantom  images.  In  the  first  approach, 
we  kept  the  contrast  between  the  first  phantom  detail  and  its  background  fixed  for  all  the 
images  displayed  in  the  study.  The  range  between  the  average  brightness  in  and  around 
the  first  detail  was  mapped  to  a  range  of  328-372  JNDs,  corresponding  to  a  range  of  73-110 
digital  driver  levels.  In  the  second  approach,  the  images  were  modified  to  have  fixed  variance 
of  the  background  noise.  Each  12-bit  pixel  value  was  multiplied  by  the  ratio  between  the 
standard  deviation  of  the  background  noise  (of  the  given  image)  and  the  standard  deviation 
of  the  background  desired  for  display.  The  new  values  of  the  average  background  and 
the  background  standard  deviation  were  selected  so  that  a  range  of  about  ±3  background 
standard  deviations  could  be  shown  by  the  8-bit  monitor  digital  driver  levels.  A  discussion 
about  the  effects  of  display  strategy  on  the  results  of  the  experiment  is  given  in  Section  IV. 
In  both  strategies,  the  linear  transformation  from  the  flat-panel  output  to  the  JND  index, 
was  image  dependent.  An  example  of  a  processed  phantom  image  is  shown  in  Figure  1. 

C.  Observer  study 

1.  Selection  and  display  of  images 

Six  medical  physicists  participated  in  the  study.  Each  observer  was  presented  with 
60  images,  consisting  of  30  for  monoscopic  viewing  and  30  for  stereoscopic  viewing.  The 
monoscopic  images  were  displayed  so  that  both  the  left  and  right  eye  saw  the  same  image, 
while  for  the  stereoscopy,  the  left  and  right  eye  images  were  different.  Each  of  the  observers 
saw  the  same  set  of  monoscopic  images  (there  were  5  images  for  each  of  6  exposure  levels). 
A  set  of  5  out  of  twenty  possible  stereo  image  pairs  for  each  mAs  was  selected  randomly  for 
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each  observer. 

The  two  sequences  of  30  images,  one  for  mono-  and  the  other  for  stereoscopic  viewing, 
were  separately  randomly  permuted  and  then  interleaved  into  a  single  sequence  of  the  form: 

...  Mn  Sn  1  'Sn-t-l  •••) 

where  Mi  and  Si  represent  mono  and  stereo  images,  respectively.  Next,  the  positions  of  the 
mono  and  stereo  images,  in  each  mono-stereo  pair,  were  randomly  changed  in  order  to  break 
the  repetition  of  mono  following  stereo,  and  to  decrease  possible  observer  bias. 

2.  Test  of  stereoscopic  acuity 

The  stereoscopic  acuity  of  the  observers  participating  in  the  study  was  tested  using  a 
standard  clinical  test,  the  RANDOT  Stereotest  (Stereo  Optical  Company,  Chicago,  IL).  The 
test  consists  of  image  pairs  with  or  without  angular  disparity.  The  images  in  a  pair  are  or¬ 
thogonally  polarized  and  the  stereo  effect  is  observed  using  appropriately  polarized  glasses. 
The  images  are  divided  into  three  groups  with  different  content  (geometric  figures,  animals, 
and  triplets  of  circles)  designed  to  assess  the  stereoscopic  acuity  in  a  wide  range  of  patients’ 
age.  The  images  in  each  of  the  groups  gradually  decrease  in  angular  disparity  and  roughly 
cover  the  range  of  20-500  seconds  of  arc  when  viewed  from  40  cm  distance.  Stereoscopic 
acuity  is  measured  by  identifying  the  image,  and  the  corresponding  angular  disparity,  be¬ 
yond  which  the  observer  cannot  distinguish  between  the  objects  with  and  without  disparity. 
Table  I  shows  the  minimum  observable  angular  disparity  for  the  six  observers. 

3.  Viewing  conditions  and  scoring 

The  C-d  experiments  were  performed  in  a  darkened  room,  with  a  black  monitor  back- 
ground.  The  distance  between  the  observer’s  eyes  and  the  monitor  was  approximately  1 
m.  The  monitor  used  in  the  study  was  equipped  for  stereoscopy  with  a  stereoscopic  gog¬ 
gles  system  (StereoGraphics,  San  Rafael,  CA).  A  transmission  box,  attached  to  the  monitor, 
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controlled  the  opening  and  closing  of  the  liquid  crystal  shutters  in  the  goggles  synchronously 
with  the  changing  of  the  image  on  the  monitor.  Therefore,  each  of  the  eyes  was  presented  the 
corresponding  image  with  a  repetition  frequency  of  56  Hz.  The  goggles  were  used  for  both 
mono-  and  stereo  experiments  to  provide  a  comparable  brightness  level  in  both  experiments. 

Several  training  sessions  for  image  scoring  were  organized  for  observers  in  order  to  es¬ 
tablish  a  uniform  decision  criterion.  The  observers  were  trained  to  inspect  the  objects  from 
largest  to  smallest  and  from  greatest  to  least  contrast.  The  objects  were  inspected  for  gen¬ 
eral  roundness  (“whether  or  not  more  than  50%  of  the  edge  was  visible”),  size  (  whether 
or  not  more  than  50%  of  the  object  was  missing”),  and  the  expected  position  in  the  detail 
array.  The  goal  was  to  prevent  misinterpretation  of  the  clustered  background  noise  as  phan¬ 
tom  details.  A  graphical  user  interface  was  developed  allowing  an  observer  to  identify  the 
smallest  detail  seen  in  a  column  by  clicking  over  the  detail  position. 

We  have  analyzed  the  relationship  between  the  inter-  and  intra-observer  variability,  by 
having  one  of  the  observers  additionally  reading  all  the  C-d  experiments  in  the  study.  Thus, 
we  could  compare  the  average  C-d  performance  of  a  group  of  observers  with  a  single  round 
of  experiments,  and  the  average  performance  of  a  single  observer  after  repeating  the  study 
multiple  times. 

For  further  comparison,  we  have  also  repeated  the  C-d  experiments  for  all  the  observers 
without  the  use  of  the  stereoscopic  goggles.  Due  to  the  112  Hz  switching  between  the  two 
stereoscopic  images  on  the  monitor,  the  observers  were,  effectively,  presented  with  an  average 
of  the  left  and  right  eye  image.  Therefore,  we  could  compare  the  performance  of  observers 
attempting  to  utilize  information  from  disparate  images  with  the  performance  of  observers 
looking  at  the  combined  images. 

III.  RESULTS 


The  observers  took,  on  average,  39  minutes  (std.  dev.  —  2.8  minutes)  to  score  a  set  of 
60  combined  mono  and  stereo  images.  Figure  2  shows  the  number  of  details  seen  mono- 
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and  stereoscopically  by  a  single  observer,  averaged  over  all  images  with  the  same  mAs.  Fig¬ 
ure  3  shows  the  difference  between  the  number  of  details  for  stereo  and  mono  images,  for 
all  the  observers.  When  averaging  over  different  observers,  we  have  included  the  results  of 
all  but  one  of  the  observers.  This  observer  reported  uncorrected  vision  in  one  eye,  corre¬ 
sponding  to  low  stereo  acuity,  and  demonstrated  a  large  value  of  the  minimum  observable 
angular  disparity  in  the  stereo  vision  RANDOT  test  (Table  I,  Observer  6).  Furthermore,  the 
difference  between  the  number  of  details  seen  stereo-  vs.  monoscopically,  for  the  observer 
in  question  (Figure  3,  points  labeled  by  triangles)  was  significantly  smaller  compared  to 
the  others.  We  have  compared  the  stereo-  and  monoscopic  performance  of  the  observers  by 
analysis  of  the  C-d  curves  and  an  analysis  of  the  SNR. 

A.  Analysis  of  C-d  curves 

For  a  low-contrast  object  with  given  thickness  t  and  diameter  d,  contrast  is  proportional 
to  t,  and  the  detectibility  is  determined  by  the  signal-to-noise  ratio,  SNR=  t  d  K\/X, 
where  X  is  the  exposure  in  mAs  and  A"  is  a  proportionality  constant.  C-d  curves  relate  the 
diameter  and  thickness  of  objects  at  the  limits  of  visibility.  Ideally,  the  ( t ,  d)  points  in  C-d 
experiments  would  form  a  line  with  negative  unit  slope  on  a  log-log  plot,  satisfying: 

logio  t  +  log10d  =  log10SNRr  —  log^K  —  -logioX,  (4) 

where  SNRr  is  the  value  of  the  SNR  at  the  threshold  of  perception. 

For  each  observer  we  calculated  their  C-d  curves  by  averaging  the  size  of  the  smallest 
visible  detail  in  each  column,  over  a  set  of  images,  e.g.,  for  the  same  exposure.  Computing 
the  average  detail  size  is  problematic  when  no  details  are  visible  in  a  given  column  for  one  or 
more  images,  since  the  contribution  of  those  images  to  the  average  is  ambiguous.  To  avoid 
this  problem,  we  chose  the  median  diameter  for  generating  C-d  curves.  Figure  4  shows  the 
average  C-d  curves  for  a  single  observer  (the  same  observer  as  in  Figure  2).  Comparison 
of  the  C-d  curves  corresponding  to  stereo-  and  monoscopic  viewing  has  been  performed 
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using  the  approximation  of  a  line  with  negative  unit  slope,  as  per  Eq.  4.  The  difference  in 
the  intercept  values  between  the  fitted  stereo  and  mono  C-d  curves  can  be  related  to  the 
effective  exposure  for  stereoradiography,  assuming  the  same  threshold  SNR  for  mono  and 
stereo  experiments,  by 

logio  du  -  logic  ds  =  logioT^  =  log10 

“s 

where  du  and  ds  are  the  diameters  of  the  smallest  visible  details  viewed  mono-  and  stereo- 
scopically,  respectively,  for  the  same  detail  thickness.  XM  is  the  exposure  used  in  acquiring 
the  projection  viewed  monoscopically,  and  Xseff/Xu  is  the  ratio  of  the  effective  stereoscopic 
exposure  (in  terms  of  object  visibility)  to  the  exposures  used  to  produce  the  individual  pro¬ 
jections.  Thus,  Xseff  is  defined  as  the  exposure  needed  to  obtain  an  image  which,  when 
viewed  monoscopically,  would  produce  the  same  C-d  characteristic  as  the  stereo  image  in 
question.  If  our  hypothesis  that  human  observers  can  combine  information  from  both  im¬ 
ages  is  true,  then  the  intercpets  of  the  corresponding  linear  fits  to  the  C-d  curves  would  be 
separated  by  y/2.  Figure  5  shows  the  linear  approximation  of  the  C-d  curves  for  a  single 
viewer. 

Figure  6  (a)  shows  the  relative  reduction  in  exposure  when  using  stereo-  rather  than 
monoscopic  viewing  normalized  with  the  hypothesized  reduction  by  half.  Plotted  are  the 
percent  values  of  Rx  -  XSe///(2  Xu)  =  [du/ (v^  ds)]2,  averaged  over  all  the  observers,  for 
both  display  conditions  for  each  per-projection  exposure  level.  If  the  hypothesis  is  true,  the 
ratio  Rx  would  be  equal  to  100%.  Figure  6  (b)  presents  Rx  for  multiple  readings  done  by  a 
single  observer  and  both  display  conditions.  For  clarity  of  presentation,  the  results  obtained 
with  different  display  conditions  are  slightly  horizontally  displaced  on  the  graph. 

B.  Analysis  of  SNR  cut-off 

The  objects  used  in  this  observer  study  were  of  low  contrast  and  imaged  on  a  detector 
with  essentially  no  correlation  in  the  noise  between  detector  elements  (a  white  noise  power 
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spectrum)15.  Therefore,  for  monoscopically  viewed  images  the  classic  Rose  model21,22  is 
applicable,  in  which  the  signal-to-noise  ratio  (SNRm)  can  be  computed  from  the  subject 
contrast  and  the  standard  deviation  in  the  pixels  values,  and  object  visibility  is  directly  re¬ 
lated  to  SNRm-  The  subscript  “M”  is  used  here  to  emphasize  that  this  is  the  SNR  applicable 
to  an  analysis  of  viewing  a  single  image.  In  Figure  7  the  number  of  objects  visible  as  a  func¬ 
tion  of  exposure  were  calculated  for  a  class  of  ideal  observers  which  would  detect  all  objects 
above  a  given  threshold  SNRxM.  Thresholds  SNRrM  of  two  through  seven  are  shown;  the 
number  of  objects  decreases  with  increasing  threshold.  Superimposed  on  Figure  7  are  the  fits 
to  the  actual  performance  of  one  observer  from  Figure  2.  For  monoscopic  viewing,  the  data 
for  this  observer  falls  between  the  theoretical  curves  corresponding  to  SNRM  thresholds  of 
5  and  6,  (in  agreement  with  the  customary  range  of  estimates  of  Rose’s  k).  For  an  observer 
capable  of  combining  information  from  stereoscopic  images  the  signal-to-noise  ratio  of  each 
object  should  increase  by  a  factor  of  y/2,  (i.e.,  SNRs  =  -s/2  SNRm),  so  that  the  object  count 
as  a  function  of  per-projection  exposure  would  correspond  to  a  monoscopic  signal-to-noise 
ratio  cut-off  reduced  by  a  factor  of  s/2-  Consistent  with  this,  the  real  observer  from  Figure  7 
does  show  an  increased  sensitivity,  as  will  be  discussed  in  detail  later. 

Following  the  Rose  model,  the  results  were  analyzed  by  identifying  a  threshold  SNR 
under  each  viewing  condition.  For  each  observer  scoring  of  the  phantom,  the  detected 
details  had  greater  SNRr  than  the  details  not  detected,  with  only  a  few  details  near  the 
threshold  not  following  this  pattern.  Thus,  for  each  observer  scoring  of  the  phantom,  the 
threshold  was  calculated  as  the  mean  of  the  smallest  SNRm  of  the  details  detected  and  the 
largest  SNRM  of  the  details  not  detected.  The  averages  of  these  per-scoring  values  over  the  5 
monoscopic  viewings  give  the  monoscopic  threshold  SNRTm  for  each  viewer  at  each  exposure 
level,  and  similarly  the  5  stereoscopic  viewings  give  SNRys .  A  hypothetical  observer  capable 
of  reducing  the  quantum  noise  in  the  phantom  images  by  using  information  from  both 
images  of  the  stereoscopic  pair  would  then  effectively  find  the  signal  to  noise  ratio  of  each 
disk  increased  by  a  factor  of  \/2,  and  thus  would  present  a  decrease  in  the  SNRm  of  the 
threshold  objects  by  a  factor  of  l/\/2.  Thus,  for  such  a  hypothetical  observer,  one  expects 
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the  ratio  RSnr  =  SNRTm/SNRTs  =  y/2.  Figure  8  (a)  shows  the  ratio  RSnr  averaged  over 
all  the  observers  at  each  exposure  level.  Figure  8  (b)  shows  the  and  the  ratio  RSnr  for 
multiple  readings  performed  by  a  single  observer.  Discussion  of  the  experimentally  observed 
results  is  presented  in  the  next  section. 

IV.  DISCUSSION 
A.  Effects  of  display  conditions 

The  transformation  from  the  flat  panel  12-bit  output  to  the  number  of  JNDs  was  imple¬ 
mented  in  an  image  dependent  fashion,  as  described  in  Section  II,  by  keeping  fixed  (a)  the 
background  contrast  of  the  largest  detail  or  (b)  the  standard  deviation  of  the  background 
noise,  in  all  the  displayed  images.  Using  the  first  approach,  images  acquired  with  high 
nominal  exposure  were  displayed  such  that  the  background  noise  range  was  limited  to  only 
a  few  8-bit  digital  driver  levels,  with  the  potential  for  quantization  distortion.  The  latter 
approach,  however,  changed  the  contrast  of  the  largest  detail  in  the  images  with  different 
exposures,  which  could  lead  to  possible  saturation  of  the  display.  Moreover,  the  percep¬ 
tual  linearization  is  a  nonlinear  transformation  that  may  affect  the  perception  of  the  noise 
patterns.  The  effects  of  these  issues  on  the  results  of  the  C-d  experiment  are  not  obvious. 

For  all  performance  measures  and  both  display  conditions,  we  observed  better  agreement 
with  theory  at  lower  exposures.  The  normalized  effective  dose  ratio  (Figure  6)  at  the  higher 
exposures,  showed  more  benefits  from  stereoscopy  for  the  images  displayed  with  fixed  noise, 
than  for  the  images  with  fixed  contrast.  The  SNR  ratio  is  near  RSnr  =  V%  at  low  exposures, 
as  expected  for  the  hypothetical  observer  combining  information  from  each  of  the  stereo  pair 
of  images.  At  higher  exposures,  the  ratio  is  less,  but  still  greater  than  one. 

The  value  of  the  normalized  effective  dose  ratio,  averaged  over  all  six  mAs  stations,  is 
J?x=89.7%±8.4%,  corresponding  to  a  reduction  in  threshold  SNR  of  i?sjvii=l-339±0.066, 
for  the  images  displayed  with  fixed  contrast,  and  Rx=  94.2%±5.0%,  corresponding  to 
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JR5ivii=1.372±0.038,  for  the  images  with  fixed  noise  (compare  with  Figure  8).  The  total 
dose  to  obtain  a  stereo  pair  D$  —  Dm /Fix  =  ^Dm/Rsnr  where  Dm  is  the  dose  for  a 
single  radiographic  projection.  Thus,  stereomammography  requires  a  dose  of  1.12±0.13  DM 
to  1.06±0.06  Dm,  for  the  respective  values  of  Rsnr- 

The  obtained  reduction  in  the  SNR  threshold  can  be  compared  with  the  results  reported 
by  Hsu  et  al.13.  In  that  study  of  stereomammographic  detection  of  simulated  objects,  the 
performance  was  evaluated  by  an  ROC  analysis.  Reported  values  of  the  area  under  the 
ROC  curve,  Az,  were  0.82  (stereo)  vs  0.74  (mono)  for  their  arrangement  experiment,  and 
0.75  (stereo)  vs  0.71  (mono)  for  their  density  experiment.  Assuming  that  Az  is  equal  to 
the  probability  of  correct  detection,  the  corresponding  reduction  in  threshold  SNR  can  be 
computed  as  the  ratio  of  stereoscopic  and  monoscopic  d!  values23,  yielding  Rsnr  =1-423  for 
the  arrangement  experiment  and  Rsnr  =1-219  for  the  density  experiment. 

We  have  performed  an  additional  validation  of  the  contrast-detail  experiment  results  by 
comparing  the  performance  measures  averaged  over  five  observers  with  the  average  perfor¬ 
mance  of  repeated  readings  by  a  single  observer  (see  Figure  6  (b)  and  8  (b)).  There  is  little 
difference  in  the  results  of  these  experiments. 

B.  Effects  of  viewing  condition 

Results  of  the  C-d  experiments  with  and  without  the  stereoscopic  goggles  are  shown  in 
Figure  9.  The  error  bars  correspond  to  the  standard  error  of  the  average  values.  For  clarity 
of  presentation,  the  2  datasets  are  slightly  displaced  on  the  graph.  Normalized  dose  ratio, 
averaged  over  all  5  observers  (solid  lines)  was  compared  with  the  performance  of  observer  6, 
with  low  stereo  acuity  (dashed  lines).  While  the  average  performance  of  the  five  observers 
with  (bold)  and  without  the  goggles  (thin)  did  not  differ  significantly,  observer  6  performed 
much  better  without  the  goggles.  This  was  expected  as  when  using  the  goggles  the  image 
shown  to  the  uncompensated  eye  of  the  low  stereo  acuity  observer  was  effectively  unused 
under  either  the  stereoscopic  or  monoscopic  viewing  conditions,  therefore  showing  basically 
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no  advantage  of  stereoscopy.  Without  the  goggles,  an  averaged  image  was  shown  to  both 
eyes,  and  the  difference  in  the  background  noise  for  the  mono  and  stereo  images  became 
visible,  without,  however,  the  advantage  of  parallax. 

C.  Effects  of  additive  detector  noise 

The  images  used  in  this  study  were  acquired  radiographically,  with  a  digital  x-ray  detec¬ 
tor.  All  images,  regardless  of  mAs,  were  x-ray  quantum  noise  limited,  rather  than  detector 
noise  limited.  Thus,  in  the  C-d  experiments  reported,  the  limiting  factor  for  detection  was 
the  x-ray  quantum  fluctuations.  Furthermore,  a  detector  entrance  exposure  of  10  mR,  typ¬ 
ical  for  the  detector  being  evaluated,  would  have  required  a  technique  of  14  mAs.  Hence, 
the  results  reported  are  likely  to  be  generally  applicable  to  digital  radiography  and  digi¬ 
tal  mammography.  Extensions  of  this  approach  can  be  envisioned  to  address  the  issue  of 
the  dose  requirements  for  tomosynthesis.  However,  in  that  instance,  careful  evaluation  of 
the  dominant  noise  source  in  the  individual  projections,  and  in  the  reconstructed  data  is 
necessary. 

D.  Other  issues  affecting  the  efficiency  of  stereoscopic  viewing 

The  following  issues  are  potential  sources  of  concern  in  our  experiment.  Since  it  was,  to 
the  best  of  our  knowledge,  the  first  C-d  study  of  the  stereoscopic  viewing  in  radiography,  no 
reference  suggestions  were  found  to  address  the  experimental  design.  Future  experiments, 
modified  to  improve  the  following  conditions,  may  give  a  more  accurate  estimate  of  the 
advantages  of  stereoradiography. 

Observers  used  the  stereoscopic  goggles  while  viewing  both  mono  and  stereo  images,  in 
order  to  keep  the  average  brightness  level  comparable  and  reduce  bias  of  having  the  mono 
images  directly  identified.  The  goggles  reduce  the  brightness  level  of  the  scene.  We  have 
not  analyzed  the  amount  of  this  brightness  reduction,  nor  its  uniformity  with  the  change  of 
brightness.  Perceptual  linearization  was  applied  without  taking  into  account  this  brightness 
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reduction,  and  this  might  be  of  concern  if  the  non-linearity  of  the  human  contrast  sensitivity 
was  not  properly  compensated.  Note,  also,  that  structured  anatomic  noise  is  not  included 
in  this  set  of  test  images,  but  its  modeling  is  considered  for  future  research.24,25  We  have 
also  begun  another  approach  to  testing  the  hypothesized  benefits  of  stereoradiography  using 
a  2-alternative  forced-choice  experiment. 

V.  CONCLUSIONS 

A  series  of  C-d  experiments  was  performed,  testing  the  hypothesis  that  the  radiation 
dose  required  for  stereoradiography  is  one  half  the  dose  for  a  single  x-ray  image  viewed 
monoscopically,  due  to  the  combination  of  images  by  the  human  visual  system.  For  images 
acquired  at  a  fixed  exposure,  more  objects  were  detected  under  stereoscopic  than  monoscopic 
viewing  conditions.  For  the  lower  range  of  the  x-ray  exposures  tested,  where  the  quantum 
fluctuations  were  easily  visualized,  the  increase  in  the  number  of  details  observed  was  the 
same  as  expected  for  a  hypothetical  observer  combining  the  two  projections  to  remove  noise. 
With  increasing  x-ray  exposure,  the  increment  between  the  number  of  objects  detected  under 
stereoscopic  and  monoscopic  conditions  decreased,  but  was  always  positive.  The  experimen¬ 
tal  results  indicate  that  a  stereoradiograph  can  be  acquired  at  a  dose  approximately  1.1 
times  that  of  a  single  projection  radiograph.  The  study  results  are  potentially  influenced 
by  the  brightness  reduction  of  stereoscopic  goggles,  background  noise  quantization,  and  the 
lack  of  anatomical  noise. 
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Fig  1  An  example  of  a  phantom  image  used  in  the  study.  The  image  was  acquired  with 
an  exposure  of  lOmAs  and  post-processed  to  compensate  for  the  effects  of  detector 
nonuniformity  and  to  adjust  for  the  non-linearity  of  human  contrast  sensitivity. 

Fig  2  Number  of  details  seen  in  the  C-d  phantom  for  two  image  display  modes:  fixed  con¬ 
trast  (thin)  and  fixed  noise  (bold).  Shown  are  the  average  number  of  details  (symbols), 
the  standard  error  (bars),  and  the  linear  fits  (lines),  for  mono-  (solid)  and  stereoscopic 
(dashed)  viewing. 

Fig  3  Comparison  between  mono-  and  stereoscopic  viewing,  (a)  Difference  in  the  number 
of  objects  seen  by  stereo-  vs.  monoscopy,  for  each  observer,  and  images  displayed  with 
fixed  contrast,  (b)  Difference  in  the  number  of  objects  seen  by  stereo-  vs.  monoscopy, 
for  each  observer,  and  images  displayed  with  fixed  noise. 

Fig  4  C-d  curves  for  a  single  observer  (observer  1),  for  monoscopic  (solid)  and  stereoscopic 
(dashed)  viewing,  and  two  display  conditions:  (a)  with  fixed  contrast  and  (b)  with 
fixed  noise. 

Fig  5  Linear  approximation  of  C-d  curves  for  a  single  observer  (observer  1),  for  monoscopic 
(solid)  and  stereoscopic  (dashed)  viewing,  and  two  display  conditions:  (a)  with  fixed 
contrast  and  (b)  with  fixed  noise. 

Fig  6  Normalized  dose  ratio,  [XM/Xs]/2  =  [(du/ds)2]/2  ,  as  a  function  of  exposure,  for 
two  display  conditions,  and  averaged:  (a)  over  5  observers  and  (b)  over  5  readings  by 
a  single  observer  (observer  3).  Bars  represent  the  ranges  of  ±  one  standard  error. 

Fig  7  Number  of  details  seen  in  the  C-d  phantom  by  an  ideal  observer,  with  various  thresh¬ 
old  SNR  values  (thin).  The  results  for  the  real  observer  from  figure  2,  the  fixed  noise 
display  and  stereoscopic  viewing  (bold  line)  are  given  for  comparison. 

Fig  8  Comparison  of  the  observers  mono-  and  stereoscopic  performance  using  the  analysis 
of  the  SNR  cut-off  values  for  two  display  conditions,  and  averaged:  (a)  over  5  observers 
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and  (b)  over  5  readings  by  a  single  observer  (observer  3).  Bars  represent  the  ranges  of 
±  one  standard  error. 

Fig  9  Comparison  of  the  normalized  dose  reduction  for  the  C-d  experiments  with  (bold 
lines)  and  without  (thin  lines)  the  stereoscopic  goggles.  Shown  are  the  results  averaged 
over  5  observers  (solid  lines)  and  for  a  single  observer  with  low  stereoscopic  acuity 
(dashed  lines). 
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Table  I  Results  of  the  RANDOT  stereo  test  for  six  observers.  Tabulated  are  the  values  of 
the  minimum  observable  angle  of  stereo  disparity  in  seconds  of  arc. 
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